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Sublimation in a Hypersonic Environment! 


SINCLAIRE M. SCALA* 


General Electric Company 


Summary 


A priori knowledge of the response of materials subjected to a 
severe aerothermal environment is essential in the space age. 
The successful design of space and re-entry vehicles demands 
that the fundamental problem of the interaction between a ma- 
terial and dissociated air be properly formulated and solved. In 
this paper, the problem of sublimation in a hypersonic environ- 
ment is considered. 

In this study of hypersonic ablation, the pertinent conserva- 
tion equations are derived and the simultaneous processes of 
diffusion, convection, and thermal exchange are analyzed for the 
vaporization of a refractory material which is subjected to the 
environmental conditions encountered during hypersonic re- 
entry. 

For simplicity, only the forward stagnation point of an axially 
symmetric body is treated. It is shown that the quantity Q*, 
called the effective heat of vaporization, which includes all heat 
absorbing or heat blocking effects, is an increasing function of 
flight speed, independent of body size, except where nonequi- 
librium vaporization effects or radiative effects appear 


Symbols 


altitude 
mass fraction of species i 
specific heat at constant pressure of species 7 


C;Cp;, frozen specific heat of the mixture 


binary diffusion coefficient 
similarity stream function 
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u/u,, dimensionless velocity 
gaseous species 


static enthalpy of species 7, including chemical 
UeGh:, static enthalpy of mixture 
t 


equilibrium constant of species 7 

frozen coefficient of thermal conductivity 
pCpD j/k, frozen Lewis Number 

liquid phase 


(2.0) , interphase mass-transfer rate 
Z w 


molecular weight of species 7 
x, M,, molecular weight of mixture 
+ 


static pressure 

omr k, frozen Prandtl Number 

Qeai./Mtot., the effective heat of vaporization 
universal gas constant 

nose radius of body 

solid state 

temperature 

x component of velocity 

stagnation-point velocity gradient 
macroscopic stream velocity 

absolute velocity of species 7 

diffusion velocity of species 7 

flight speed 

mole fraction of species 7 

body-oriented coordinate system 

pure element 

m, /Mtot., gasification ratio 

similarity variables 

T/T. 


coefficient of viscosity 


., dimensionless temperature 
vaporization coefficient 
emissivity 

Stefan-Boltzmann constant 
density of gas 


Subscripts 


calorimeter 

outer edge of the boundary layer 
equilibrium 

gas cap (shock layer) 

ith species 

vaporizing species 

liquid phase 

stagnation point 
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tot. = total 
w = wall, interface 
n = denotes partial differentiation with respect to 7 
0 = sea level, standard 
oo = upstream of shock 


Introduction 


; 
A BODY immersed in a high-temperature, high-ve- 

locity stream behaves initially as a pure heat sink. 
When the surface reaches a critical temperature, selected 
particles at favorable sites will become detached from 
the surface lattice structure. And then, during the 
microscopic collision processes which follow, a number 
of particles are readsorbed by the surface, while others 
are transported away by convection and diffusion in the 
fluid stream. 

Sublimation is the name given to the process of 
vaporization directly from a solid phase into the gas 
phase. In general, however, many materials exhibit 
an intermediate molten state or liquid phase when sub- 
jected to a high heat flux, even at ordinary pressure 
levels (bp ~ 1 atm.). And, depending on the magnitude 
of aerodynamic and body forces which act at the surface, 
the viscous liquid layer which forms will either remain in 
place or be swept downstream. Thus, if one introduces 
the concept of a gasification ratio, defined by 

mass rate of vaporization 


(1) 


local total rate of mass loss 


it is clear that sublimation corresponds to the limiting 
case of ! approaching unity. 

In the past, the theory of mass transfer into the high- 
speed boundary layer!~* has been confined mainly to the 
study of single component fluids. More recently, 
studies of the effect of the injection of foreign gases into 
the boundary layer have also begun to appear.‘~’ A 
current problem arising from the technology of re- 
entry from space is the analysis of mass transfer due 
to chemical reactions, and this has focused attention on 
the study of heterogeneous chemical reactions such as 
pyrolysis, depolymerization, combustion, and vaporiza- 
tion®—'4 for the conditions of hypersonic heat transfer. 

In order to confine the region of thermal degradation 
and structural failure to a narrow zone near the surface, 
it is essential that the ablating material have a low co- 
efficient of thermal conductivity, a high viscosity in the 
molten state (if one forms), and a substantial heat of 
phase change. 

From both a practical and theoretical point of view, 
it is interesting to consider, as a class, quartz-like ablat- 
ing refractories, which not only meet the requirements 
demanded of ablating materials, but also present an 
interesting problem in multicomponent fluid flow. 

Of particular interest is the sublimation of a refrac- 
tory dioxide, written symbolically ZO.(s), where the 
symbol Z denotes the pure element. If vaporization 
occurs, then the foreign gaseous species which vaporize 
from the surface include the element oxygen, the mon- 
oxide ZO and the dioxide ZO», where it is noted that one 
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half mole of oxygen molecules accompanies each mole 
of ZO. That is, the surface vaporization reactions can 


be written 


ZOA(s) = ZO(l) = ZO(z) + 1/2 Ow (2) 
ZO2(s) = ZO2(l) = ZO2(g) (3) 


In words, these equations state that the original solid 
dioxide material softens and forms a quasi-liquid molten 
layer which vaporizes, producing an injected gas which 
contains one of the components of dissociated air. Since 
the mass-transfer rate is a strong function of the com- 
position of the gas at the surface, the equilibria of such 
a system must be carefully considered.” It will be 
shown that when the heterogeneous reaction kinetics 
are assumed to be infinitely fast, so that the gas at the 
surface is maintained in local thermochemical equi- 
librium—i.e., the overall rate process is assumed to be 
then the equilibrium constants’® 
‘equilibrium rate of 


diffusion controlled 
of the oxygen-nitrogen system, the 
vaporization” and the fact that the surface is imperme- 
able to the element nitrogen, all taken together, provide 
a sufficient number of constraints at the surface to 


‘ 


make the gaseous composition determinate. Note that 
once the composition of the gaseous mixture at the inter- 
face is known, all other physicochemical quantities 
follow by quadrature. 

In analyzing the problem of melting and vaporization 
it has been convenient in the past to assume the ex- 
istence of a characteristic surface ‘“‘vaporization tem- 
perature’ which is taken as a property of the material 
alone. In this more general study, the surface tempera- 
a priori,’ but rather 


sé 


ture is not assumed to be known 
T, is obtained by applying the conservation laws for 
the interphase mass and energy transfer at the interface 
between the gaseous and condensed phases, which 
yields the solution to the problem. 


Rate of Vaporization into Gas Phase 


Scala and Vidale recently made a detailed study of 
vaporization processes in the hypersonic laminar 
boundary layer‘ in which it was assumed that gas 
phase reactions were frozen. It was found that, for all 
environmental conditions, the interphase mass transfer 
mi, Was directly proportional to the mass fraction of the 
vaporizing species evaluated within one mean free path 
of the surface C;,,, and that the proportionality con- 
stant depended on the environmental conditions. The 
functional form of this relationship was given as 


tity V Ra! Crm = fi(Alt., Vor Tr) (4) 
Numerically, this was expressed as 


lite Re/ Cry = [1.233 — 6.6667 X 10-57,,(°R.)] X 


— [4.32 — 3.97 6 
107 | 3.97 X 10-6 Alt fl) x 
(Vv [0.9976 — 1.201 x 10-6 Alt. (ft 11] (Ib. /ft.3/2 sec.) (3) 
2 , (Ib. /ft.*/? sec. : 
at the forward stagnation point of an axisymmetric 


body, where I’. is in ft. /sec. 
Thus, a knowledge of the driving force C;.,, immedi- 
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SUBLIMATION IN A 


ately yields the rate of mass transfer m,, provided it 
can be shown that the frozen gas assumption made in 
reference 14 is compatible with the problem at hand. 
Some discussion is therefore required concerning the 
magnitude of the rates of homogeneous and hetero- 
geneous chemical reaction. 

In solving the gas phase boundary-layer equations, 
it has usually been assumed, for convenience, that 
either (a) homogeneous chemical reactions are every- 
where frozen, or (b) the gas is everywhere in “‘local 
thermochemical equilibrium.’ The first of these as- 
sumptions implies that the rates of chemical reaction 
are everywhere zero in the gas phase. The second as- 
sumption implies that the gas composition is deter- 
mined by the simultaneous solution of a system of 
equations representing the equilibrium constants of the 
governing reactions. While not always realistic, these 
assumptions are nevertheless quite useful in that these 
chemical constraints yield the two limiting possibilities 
in the realm of nonequilibrium boundary-layer flow, 
which presumably bracket the actual solution. 

In a study of the nonequilibrium dissociated binary 
mixture air laminar boundary layer, the author ob- 
tained the anticipated general result that, for a given 
level of dissociation in the free stream, the concentration 
of atomic species at a surface having finite catalytic 
efficiency depends on both the homogeneous and 
heterogeneous rates of atom recombination.” Further- 
more, regardless of the specific combination of these gas 
phase and surface reaction rates, for given free-stream 
conditions, the heat transfer to the surface varies in- 
versely with the mass fraction of air atoms in the gas at 
the surface. This result was obtained because the heat 
transfer to a surface depends primarily on the enthalpy 
difference, including chemical, across the boundary 
layer, and hence for given free-stream conditions de- 
pends, to the first order, only on the precise composi- 
tion of the gas at the surface. Thus, for example, if the 
gas is in equilibrium with the surface, it is generally 
immaterial as to whether this equilibrium composition 
is obtained as a result of the magnitude of the homo- 
geneous or of the heterogeneous reaction rates 

Now, a fully catalytic surface is defined as one which, 
if placed in contact with a gas, will produce an equi- 
librium composition. That is, the gaseous particles in 
the first few mean free paths which continually undergo 
collisions with the surface, have sufficient dwell time 
so that the gas rapidly approaches an equilibrium 
composition. This concept can also be stated as 
follows: the gas in contact with a fully catalytic surface 
reacts chemically at virtually infinite rates, such that 
the composition is in “‘local equilibrium” within the 
first mean free path from the surface. Thus, by as- 
suming that a surface is fully catalytic, which is 
equivalent to assuming that the heterogeneous reaction 
rates are virtually infinite, one can compensate for the 
assumption of zero homogeneous reaction rates, since 
the heat transfer should then be the same as for the 
assumption of “local equilibrium’ throughout the 


boundary layer. 


HYPE 


RSONIC ENVIRONMENT 3 
In extending this concept to the case of vaporization, 
it is first important to note that when the interphase 
mass transfer is finite, the mass fraction of vaporizing 
species at the surface cannot be identically equal to the 
equilibrium mass fraction of vaporizing species. This 
follows because under ‘“‘true’’ equilibrium conditions 
the condensation rate equals the vaporization rate, and 
the net rate of mass transfer is zero. That is, kinetic 
theory’ yields, as the net interphase mass transfer, 


tity = (aM ypP/V24MRT z)(Cieg — Cdw (6) 


This relation must be compatible with Eq. (4). As dis- 
cussed in reference 14, when the vaporization coefficient 
a is sufficiently large, the resistance of the material to 
the vaporization process is small compared with the fluid 
dynamic resistance of the diffusion and convection 
process. Thus, when a 2 aerit., the overall rate deter- 
mining step is the diffusion process and then C;..,,, = 
Cx, in which case one may utilize the equilibrium con- 
stant of the vaporization reaction in determining the 
composition of the gas at the interface. 

Thus, when we make the assumption that the overall 
process is diffusion controlled (large a) in this study, 
this is also equivalent to the assumption that the surface 
is fully catalytic and that the gas is in local thermo- 
chemical equilibrium. This, in turn, implies that the 
rates of chemical reaction in the gas phase will play a 
minor role in the overall exchange of mass and energy 
in the boundary layer, in which case Eq. (5) applies, 
provided only that the vaporizing species have transport 
properties (collision diameters, force constants, and 
molecular weights) which are not strongly dissimilar 
from air molecules'® and, hence, result in Lewis Num- 
bers of approximately unity. 


Conservation of Mass, Momentum, and Energy 
at the Interface 


The conservation of mass during the ablation process 
requires that the total mass ablated from the solid 
phase must appear either in the liquid phase or as 
gaseous products of vaporization which enter the gas 
phase at the interface between the gas and liquid phase 
boundary layers. 

The conservation of momentum requires that the 
aerodynamic shear stresses generated by the viscous gas 
phase boundary layer, which flows over the surface of 
the material under the influence of an external pressure 
gradient, be balanced by the strain of the condensed 
If the condensed phase is a solid, its elastic 
If the condensed 


phase. 
modulus will determine the strain. 
phase is in a molten (liquid) state, the deformation of 
the surface fluid elements results in viscous flow. 

The conservation of energy dictates that the aero- 
dynamic heat transfer generated by the viscous interac- 
tion between the stream of dissociated air and the sur- 
face of the material be either blocked or absorbed during 
the degradation of the surface layer of the material. 

Although a space vehicle encounters time-dependent 
during re-entry, in this 


environmental conditions 
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Fic. 1. Coordinate system and profiles for vaporizing boundary 
layer. 


analysis we will treat the quasi-steady ablation prob- 
lem, in which the free-stream conditions are held fixed 
and, hence, partial derivatives with respect to time are 
zero. Thus, an observer stationed at the interface will 
see a steady flow of gas m, from the condensed phase 
into the gas phase boundary layer (see Fig. 1). For the 
sake of generality it is assumed that a liquid layer exists. 
Thus, when 0 < I < 1.0, a stream of molten material 
flows away from the stagnation point at a steady rate ?2;. 
For simplicity, but with no loss of generality, it will also 
be assumed that the molten layer has the same thermo- 
dynamic properties as the solid phase. That is, the 
liquid layer is merely a softened solid layer which can 
flow under the influence of external forces, and hence no 
sharp phase change is cbserved as the solid passes into 
the molten state. A mass balance at the interface (see 
Fig. 2) requires that 


Mot. = My + Mm, = Trmrtor. + (1 —T)meot. (7) 


where 7;,¢. is the total rate of mass loss of solid material, 
Mi» is the interphase mass transfer, and m, represents the 
convection of mass in the liquid layer. (Note that it has 
already been assumed that the liquid layer thickness 
is not a function of time—i.e., liquid does not accumu- 
late at the stagnation point.) 

In considering the absolute heat transfer at a chemi- 
cally reacting interface it can be shown that for the 
gaseous layer at a stagnation point, the heat transfer 
into the surface is given by” 


[-—Qw] = 
[R(OT/Oy) — YE piVihi — mh + yu(du/dy) + Qraa.],p 
(8S) 


which includes the effects of thermal conduction, dif- 
fusion, convection, shear work, and radiation. At a 
stagnation point the shear work vanishes, in the absence 
of slip flow. 

The radiation term which can be written in the form 


Cas 


oles] s* — enol »') (9) 


*w 


represents the net radiation at the surface, due to the 
radiation toward the body from the gas cap and the re- 
radiation from the surface. 

A quasi-steady-state energy balance written for a 
control volume attached to the interface, and passing 
through the gaseous, liquid, and solid phases (see Fig. 2) 


requires that 
[-—O,] + Mot A(S) — mh, = U0 (10) 


At this point it is convenient to introduce an energy 
transfer function (Q),,, which is defined as the sum of the 
gas phase conduction and diffusion terms—.e., 


On — [R(OT Oy) re » a piV hile (11) 


so that at a stagnation point the absolute heat transfer, 
Eq. (8), becomes, in the absence of slip flow, 


[—Qe] — Ow = Mh 7 Crna. (12) 


In this expression, h, is the enthalpy of the gaseous 
mixture, including free-stream components which dif- 
fuse to the wall, products of vaporization, and products 
of chemical reaction, evaluated at the surface tempera- 
ture 7,,. (It should be noted that Cn. as defined above, 
is identical with the absolute heat transfer [—Q,,] 
when mass transfer or radiation effects are absent.) 
Upon postulating that the alteration in the function Q,, 
due to mass transfer is directly proportional to the mass- 
transfer rate, even for large rates of vaporization, one 
can write the Taylor expansion (holding the surface 


temperature constant), 


Ow = (Qu) sitgeg — (AQ/ Ami) wtiry (13) 


since 


(Que) sitapno (14) 


(Ou ) my—0 = 


Upon introducing the symbol Q..:. to represent the 
net heat transfer to a non-ablating calorimeter having 
the same surface temperature, catalytic efficiency, and 
surface emissivity as the ablating surface, one can write 


eat. = (Qu) mrp + gle,t 5" _ 6.4") (15) 


9 


Then, upon introducing Eqs. (7), (8), (11), (12), (13), 
and (15) into Eq. (10) and dividing by mtot., one im- 
mediately obtains the following expression for the effec- 
tive heat of vaporization: 

Q* = Oni: Mot. 


; T'}(AO/Anit) wo + [tw — h(s)] + (16) 
(1 — P)/T] la — h(s)]}) 


I 


In the above equation, I is the gasification ratio, which 
has been previously defined as the fractional part of the 
total mass loss which vaporizes, /1(s) is the enthalpy of 
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the cold solid ZO2(s), evaluated within the interior of 
the material, and /, is the enthalpy of the liquid ZO,2(/) 
layer evaluated at the velocity weighted average tem- 
perature of the liquid layer. If the liquid layer is very 
thin, it is reasonable to assume that the enthalpy /,,, 
evaluated at the interface temperature is substantially 
the same as the enthalpy of the liquid evaluated at its 
velocity weighted mean temperature. This form shows 
clearly that the magnitude of Q* depends upon the 
contributions of three coupled effects—i.e., 


T'(AQ Am), is the heat blocking action term due to 
boundary-layer thickening 


['[h, — h(s)] is the heat absorbing action term due to 
phase changes and gas heating upon injection 


(1 — I)[h, — h(s) ] is the convection of heat in the liquid 


phase 


and also shows the strong dependence of Q* upon the 
gasification ratio I. 

In a study’ *! of the effects of the transfer of air 
molecules into the dissociated hypersonic laminar 
boundary layer, the author found that the ‘‘effective- 
ness quotient” (AO Am),, of the mass-transfer process 
was given by the following empirical correlation 


formula: 


(AQ/ Anit)» & (30/M x) (240 + 0.48(h, — h)), 
B.t.u./Ib. (17) 


when /), the enthalpy, is expressed in B.t.u./lb. and J, 
is the mean molecular weight of the gaseous mixture 
evaluated at the surface of the condensed phase. Since 
the stagnation enthalpy /, can be of the order of 12,000 
B.t.u./lb. during hypersonic re-entry, it is clear that the 
heat blocking term '(AQ/Am). lends a significant con- 
tribution to Q*. Although not precise, we will assume 
here that the heat blocking action due to the mass trans- 
fer of the vaporizing species is given closely by Eq. (17), 
which as noted had been developed rigorously for air 
molecule mass transfer. 

The term I'{h, — h(s)] is related to the energy re- 
quired to break the intermolecular bonds during the 
vaporization process, and can be of the order of 8,000 
B.t.u./lb. for some refractory oxides, if pure sublima- 
tion occurs. 

The term (1 — I')[h,; — h(s)] represents the heat re- 
quired to bring each unit mass of solid to the molten 
state. 

When T = 0, pure melting occurs and Q* becomes 
essentially the heat of melting (in the steady state). 
Since [h,; — h(s)] is very small compared to the other 
terms, it is clear that pure melting represents an ineffi- 
cient utilization of the material. 

Examination of the right-hand side of Eq. (16) indi- 
cates that, for given environmental conditions, the 
unknowns in the problem are the gasification ratio I, 
the surface temperature 7,,, and the gas composition at 
the interface. A knowledge of these three quantities 
enables one to completely evaluate the enthalpies and 
the effectiveness quotient. Although sublimation 

















[= Mw / ry rot 


Fic. 2. Control volume for energy balance at stagnation point 


corresponds identically to the case Tl = 1.0, we have 
indicated that the entire range 0 < I < 1.0 will be con- 
sidered. To bring out the essential features of the re- 
sponse of a vaporizing material to hypersonic heating, T 
will be treated as an independent variable whose magni- 
tude may be established a postericri by matching gas 
and liquid phase boundary-layer behavior as was done 
by Scala and Sutton.'* Thus, for discussion purposes 
the unknowns are 7°, and the gas composition at the 
surface. 

Let us now consider the left-hand side of Eq. (16). 
For given environmental conditions, the numerator 
depends only on the as yet unknown surface tempera- 
ture, since the heat transfer at the stagnation point of 
a calorimeter can be written symbolically in the 


form? 23 


Qca.WR p = fo(Alt., v. ’ ee €u VR B ) (18) 


If the mass-transfer process is diffusion controlled, 
then we may introduce Eq. (4) rather than Eq. (6) into 
the denominator on the left-hand side of Eq. (16), which 
can then be written in the form 


rittot.V R p ™ tity R B ‘T= 


(Ci./T)-f,(Alt., V., Ty) (19) 


This indicates that the unknowns are the effective mass 
fraction for vaporization C;,, (i.e., the gas composition), 
and the surface temperature 7, (at a given value of I). 
Upon taking the ratio of Eq. (18) to Eq. (19) one ob- 
tains 


Qea rf.(Alt., V., Tw, tuV Ry ) 


ae = 
¢ Tx) 


= (20) 
Mot. Crp fi(Alt., J 


Thus, it is clear that if the mass-transfer process is 
diffusion controlled rather than rate controlled, such 
that mi, is independent of a, and if radiation effects do 
not appear, then the effective heat of vaporization at the 


stagnation point is independent of the size of the 


body—i.e., the nose radius R , does not appear either 
explicitly or implicitly. (Note that Eq. (20) shows that, 
for given flight conditions, Q* may be expected to be a 
linearly increasing function of I.) 

In the following section we will discuss the physico- 
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Fic. 3. Equilibrium compositions during vaporization. 


chemical constraints which lead to the determination of 
the gas composition and the surface temperature. 


Physicochemical Constraints 


If the composition of the gas at the interface can be 
related to the surface temperature, then Eq. (16) be- 
comes effectively an equation in a single unknown 
namely, 7, (at a given value of I). 

We will begin with the assumption that the pre- 
dominant vaporizing species are the monoxide ZO and 
the element oxygen—i.e., it is assumed that the 
vaporization reaction is given by Eq. (2) and Eq. (3) is 
discarded. Thus the total number of species present in 
the gas phase boundary layer is taken as five, and in- 
cludes the major products of dissociated air, namely, O, 
Oo, N, No—1.e., nitric oxide is neglected—and ZO. 

The equilibrium constant for the vaporization proc- 
ess corresponding to Eq. (2) is given by 


K pzo = (Pz) (Po,) 2 = (Xz0) (X,,) — (21 ) 


where Pzo is the partial pressure of ZO and P,, is 
the partial pressure of oxygen molecules, at the surface. 
It is remarked that X,, includes not only the con- 
tribution of the vaporizing oxygen molecules but the 
oxygen already present in the dissociated air boundary 
layer. It is clear then that the presence of free oxygen 
in the boundary layer actually inhibits the vaporization 
process, since the magnitude of X,, is reduced as Xo, 
increases, at a given value of the interface temperature. 

The equilibria of the pure oxygen and pure nitrogen 
system is given by the dissociation reactions 
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Os = 20 (22) 
Nz = 2N (23) 


for which the equilibrium constants are 
K pg = (P,)° Fa = HA) Xo |P 
Kpy = (Py)?/Py, = ((Xy)?/Xx,]P 


(24) 
(25) 


Thus, Eqs. (21), (24), and (25) may be utilized in de- 
termining the composition of the gas which is assumed 
to be in equilibrium with the surface. However, at 
moderate pressures and at the maximum surface tem- 
peratures of interest; ie., 7, ~ 5,000°R.; Kp, is 
small and hence Xy, = 0. That is, there are virtually 
no nitrogen atoms present at the surface, consequently 
Eq. (25) may be eliminated. At this point, one may 
conclude that the gas at the wall consists of a mixture 
of O, Oo, No, and ZO, while the free stream components 
are O, Oo, N, and No. 

Upon introducing Dalton’s law 

P=)>P;; > X;= 1.0 
i 


? 


(26) 


Eqs. (21), (24), and (26) may be combined to yield the 


single algebraic equation in Xo,,, 


(Xo.,,)° 2 + (Xo,,,) (Keo ry} 2 + (K pzo P3/2) = 
(Xo,,,)/2(1 — Xx,,) (27) 


“ul 


However, this equation, taken alone, does not permit 
the determination of the composition of the gas at the 
wall, since a priori knowledge of the mole fraction of 
nitrogen molecules at the surface would be required. 
Eq. (27) isa single equation in two unknowns, Xo,,, and 
X Nop. 
centrations could be determined immediately by first 
solving Eq. (27) and then using back substitution into 
Eqs. (21) and (24). It is noted parenthetically that one 
is not permitted to assume, for example, that the con- 
centration of the element nitrogen is constant across the 
boundary layer. For that matter, unlike flow in a nozzle, 
the ratio of the mass fraction of the element oxygen to 
the element nitrogen need not be conserved throughout 
the boundary layer, which is a diffusion field, since 
preferential diffusion of the various species can and does 


If XN2, were known, the remaining three con- 


occur. 

A further complication arises from the fact that the 
simultaneous solution of Eqs. (21), (24), and (27) is 
double-valued. Thus, for example, at a total pressure 
of 1.79 atm., and at four different temperatures in the 
vicinity of 5,000°R., one obtains the results shown in 
Fig. 3. For each value of Xwn»,, taken a1bitrarily in 
the range 0 < XwNo,, < 1.0, there exist two sets of 
corresponding Xo, X0:,, and Xzo,,. Inspection of 
Eqs. (21), (24), and (27) indicates that Xo, and X0»,, 
increase together, while Xzo,, decreases if either Xo,, 
or X0»,, increases. Hence, the upper portions of the 
Xo, and X0»,, curves taken together with the lower 
portion of the Xzo,, curve constitute the first set, while 
the lower portions of the Xo,, and X02, curves taken 
together with the upper portion of the Xzo,, curve 
constitute the second solution. 
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SUBLIMATION IN A HYPERSONIC ENVIRONMENT 


Although the second solution can be discarded by 
means of a physical argument regarding the magnitude 
of the mass-transfer rate as the surface temperature 
drops, the equilibrium composition is nevertheless in- 
determinate until one additional consideration is intro- 
duced. This additional equation may be derived by 
noting that the vaporizing surface is impermeable to the 
element nitrogen. 

The mass transfer of vaporizing species into the 
boundary layer drives the free-stream gaseous com- 
ponents away from the wall, and, hence, the composi- 
tion of the gas at the surface is altered from that of free- 
stream dissociated air. The mass fraction of the ele- 
ment nitrogen at the surface is thus considerably smaller 
than in the free stream, and this difference in the con- 
centration of nitrogen across the boundary layer causes 
nitrogen molecules to diffuse toward the surface. 

Since no nitrogen can penetrate the surface, the mass 
flux of the element nitrogen which is carried away from 
the surface by the mean mass motion must be equal and 
opposite to the diffusion flux of nitrogen toward the sur- 
face. Thus, the condition of the impermeability of the 
surface to the element nitrogen leads to a relation be- 
tween the interphase mass transfer and the mass frac- 
tion of the element nitrogen at the wall Cy»,, [see Ap- 
pendix, Eq. (A-7)], which can be written symbolically : 


tity Rp = fa(Alt., V., Tey C2,.) (28) 


Now the net mass transfer due to vaporization into 
the boundary layer depends on the contribution of the 
monoxide ZO and the element oxygen—that is, 


My = (Mzo + Mo, + Mo)w (29) 


Examination of Eq. (2) indicates that stoichiometry re- 
quires that the interphase mass transfer 7, be related 
to the mass flux of the monoxide by 


My = [1 + (Mo,/2M zo) \m zo, (30) 


Since the mole fraction is related to the mass fraction 
by 

C, = X,M,/M (31) 

Eq. (30) implies that the effective mass fraction for mass 

transfer, to be utilized in Eq. (5), is related to the mole 


fraction X zo, by means of 
Crp = (1 + (Mo,/2M zo) (Xz0,Mz0/ Mu) (32) 
where 
M, = (16Xo + 32Xo, + 28Xy, + MzoXzo)w (33) 


If the mole fraction of nitrogen at the wall is known, 
then eqs. (27), (21), and (24) yield the composition at 
the surface. When these values of the mole fractions 
are substituted into Eqs. (32) and (33) one obtains C;,,,, 
which may be utilized in Eq. (5) to determine the inter- 
phase mass transfer. Clearly then, Eqs. (5), (32), and 
(33) may be considered equivalent to the functional re- 
lationship 
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tity Rp = f(Alt., V., Te, Cre (34) 


-u 


Thus, for given flight conditions, and at a given surface 
temperature, 7), the unknown mass fraction Cn», can 
be eliminated between Eqs. (28) and (34). That is, the 
gas composition at the wall is now a known function 
of the surface temperature 7°, since the nitrogen com- 
patibility relation gives us a means of determining 
CN>,, at each value of the surface temperature, and the 
use of the equilibrium constants yields the remaining 
composition. 

In Fig. 4, we have shown graphically how the mass 
fraction CN»,, is established at each value of 7,,. The 
curve which descends from left to right satisfies Eq. 
(A-7) and it is noted that in the range of surface tem- 
peratures of interest only a single curve is obtained. 
The shorter curves which ascend from left to right 
satisfy Eqs. (5), (32), and (33), and were evaluated 
for four values of the surface temperature in the range, 
4,800 < 7, < 5,200°R. The intersections of the two 
families of curves uniquely relate the gas composition 
to the surface temperature. 

Thus, having related the gas composition to the sur- 
face temperature by means of the simultaneous solution 
of the equations for (a) the equilibrium constants, (b) 
Dalton’s law, (c) the interphase mass transfer, and (d) 
the impermeability of the surface to the element nitro- 
gen, one may return to the solution of Eq. (16), which 
may now be considered as an equation in the single 
unknown 7’,,. 

Note that in a rigorous treatment of vaporization it is 
not correct to assume the existence of a constant abla- 
tion temperature, since the determination of the varia- 
tion of the surface temperature with flight conditions 
yields the solution to the problem. 


Transport and Thermodynamic Properties 


For the numerical calculations, quartz was selected 
as being typical of the class of materials denoted by the 
symbol ZO. Thus, the molecular weights were taken 
as Mz = 28 and Mzp = 44. 

During the evaluation of Eq. (A-7), it was assumed 











that the Prandtl Number of the mixture and the Lewis 
Number of the diffusing molecules could be evaluated 
by treating the gas as a binary mixture of air atoms and 
air molecules;'® i.e., treating the contribution of ZO as 
equivalent to air molecules. The frozen specific heat 
Cp,, was evaluated in terms of the contributions of air 
atoms and air molecules,!? where 


C, = LD CiCy and C,,5 = 0.21 B.t.u./Ib.°R. 
i 


Upon utilizing a heat of formation” (based on a ref- 
erence temperature of 537°R.) for the monoxide of 


Ah’... = — 1,147 B.t.u./Ib. (35) 


as well as the data of references 16 and 19, the enthalpies 
of the four gaseous species present at the surface were 
expressed as a linear function of the surface temperature 
as follows: 


ho = 6,493 + 0.31057,,(°R.), B.t.u./Ib. (36) 


ho, = —200 + 0.287,,(°R.), B.t.u./Ib. (37) 
hy, = —280 + 0.327,,(¢R.). B.t.u./Ib. (38) 
hzo = —1,260 + 0.217,,(°R.), B.t.u./Ib. (39) 


The heat of formation of both the solid and molten 
quartz, based on a reference temperature of 537°R., was 
taken" as 


Ah’ = —6,200 B.t.u./Ib. (40) 


JSiO2 


in which case the enthalpy of the condensed phase may 
be written 
hzo,(s) — hzo,(l) = — 6,350 + OS77T R.). 
B.t.u./Ib. (41) 
The equilibrium constants were obtained by curve 
fitting thermodynamic data provided by Vidale’ and 
Linevsky.'® Values of the equilibrium constants which 
correspond to Eqs. (21), (24), and (25), respectively, 
are the following: 


1.020 % 10° 


| 10 K an = 13.90 <— : (42) 
Og) PZ ) TOR.) 
4.68 X 104 
logip Keg = 6.80 — — (43) 
—e T(°R.) 
9.00 * 10% 
login Key = 7.00 — (44) 
— FOR.) 


Thus, at any value of the temperature and pressure, 
Eqs. (42) and (43) can be utilized in the determination 
of the mole fractions X zo,, Xo,, and X0»,, and then 
the enthalpy of the gaseous mixture at the surface 
follows directly from the definition 


X.M, ; 
Res = (x Cih eas bP — 7] i | (45) 


1 


Dependence of Effective Heat of Vaporization 
Upon Environmental Conditions 


In the preceding sections equations were presented 
for (a) the conservation of mass, (b) the conservation 
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of energy, (c) the interphase heat transfer, (d) the inter- 
phase mass transfer, and (e) the heat blocking action of 
the mass transfer process, and these were combined 
into Eq. (16). 

In addition, Dalton’s law was combined with the 
equations for the equilibrium ZON system and with the 
nitrogen compatibility equation, in order to relate the 
gas composition to the surface temperature. 

Finally, a sufficient number of auxiliary relations 
were introduced for the property data, including (a) the 
transport properties, e.g., the Prandtl and Lewis Num- 
bers; (b) the thermodynamic properties including 
heats of formation, specific heats and enthalpies; and 
(c) the equilibrium constants, so that all of the terms 
in Eq. (16) could be evaluated as a function of 7, at 
any flight speed and altitude in the hypersonic regime. 

In this section, we will therefore proceed with the de- 
termination of Q*. The latter is a most useful quantity, 
since once the effective heat of vaporization of a given 
material has been determined as a function of environ- 
mental conditions, the aerothermodynamic designer 
can simply predict the total mass loss for a particular 
flight speed and altitude by dividing the calorimeter 
heat-transfer rate (to an equivalent surface having zero 
mass transfer) by Q*. It is noted here that for corre- 
lation purposes, flight speed and altitude effects are re- 
flected through the separate effects of stagnation en- 
thalpy and stagnation pressure. 

For prescribed environmental conditions, and a given 
value of I, the left-hand side of Eq. (16) is a strongly 
decreasing function of the surface temperature. This 
follows because the numerator, which is the heat trans- 
fer?’ to a calorimeter, decreases as the surface tempera- 
ture rises, and further, the denominator, which is the 
ablation rate, is an exponentially increasing function of 
surface temperature, for constant T. Now, for pre- 
scribed environmental conditions, the right-hand side of 
Eq. (16), which is a strong function of I’, is relatively in- 
sensitive as the surface temperature. Thus, Eq. (16) is 
satisfied uniquely at only one value of the surface tem- 
perature, if lis held constant. Then, back substitution 
yields the interphase mass transfer m,,, the ablation rate 
Nitto. and the effective heat of vaporization, Q*. A typical 
solution, at a flight speed of 18,000 ft./sec. and an altitude 
of 120,000 ft., is given in Fig. 5, in which 7,,, tityW R By 
Tittor. WR pg, and Q* are all shown as a function of I 
(assuming €,, = 0). 

It is immediately remarked that due to the exponen- 
tial dependence of m, on T,,, the surface temperature 
is relatively insensitive to the value of [. In fact, it is 
seen that, for a vaporizing material with a low thermal 
diffusivity, even when T is low (i.e., T ~ 0), the surface 
temperature can rise rapidly to a maximum value which 
is only about 200 or 300°R. below that which is ob- 
tained for pure sublimation (Tf = 1.0), without result- 
ing in any interphase mass transfer. That is, when the 
molten layer runs off without vaporizing, the surface 
temperature can nevertheless be quite high. 

On the other hand, the total mass loss is a very strong 
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Fic. 5. Dependence of mass transfer, effective heat of vaporiza- 


tion, and interface temperature on T. 


function of I’, and it was found that the effective heat of 
When T = 
1.0, the total mass loss is a minimum, and is, in fact, 


vaporization Q* increases linearly with I. 
equal to the interphase mass transfer ,._ This implies 
that the heat blocking action due to the thickening of 
the boundary layer and the heat absorbing action due 
to the phase change during vaporization and subse- 
quent heating of the gas within the boundary layer are a 
considerably more efficient means of reducing the heat 
transfer to the interior of the solid than convection in 
the molten layer of materials which do not have a sig- 
nificant heat of phase change at melting. 

When TL is identically zero, Q* is a minimum and 
convection in the liquid layer merely absorbs the heat 
necessary to bring the solid up to the molten state. 


This is shown by setting T = 0 in Eq. (16)—i.e, 
lim Q* = h, — h(s) (46) 
r—0o 


where /;; may be evaluated at the surface temperature, 
and h(s), the enthalpy of solid ZOs, is evaluated at the 
interior temperature of the cold solid. Let us now re- 
turn to the right-hand side of Eq. (16). 

For constant I, the first of the three terms increases 
almost linearly with the square of the flight speed. 

The second term is, in general, dominated by the 
The magnitude of h(s) is 
a constant which is evaluated at a low temperature 


enthalpy of the solid, h(s). 


(540°R.) and is very nearly equal to the heat of forma- 
tion of solid ZO, (—6,200 B.t.u./Ib.). 

The enthalpy of the gaseous mixture at the interface 
h, was evaluated by taking the summation: h, = 
(> Cihi)w. When T = 0, hy is equal to the enthalpy 


of equilibrium air evaluated at the interface tempera- 
ture, since the mass fraction of the injected species is 
zero for pure melting. WhenT = 1, 7), is a maximum, 
C zc, is also a maximum and consequently, /, is a 
minimum, since / zo is negative relative to the com- 
ponents of dissociated air. That is, at constant stagna- 
tion pressure the increased concentration of ZO, which 
has a large negative heat of formation, more than com- 
pensates for the slight rise in surface temperature as T 


varies from zero to unity. However, since there is a 
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pronounced dependence of surface temperature upon 
stagnation pressure, h,, correlates favorably as a linearly 
increasing function of the stagnation pressure. In this 
case, the increased contribution of the positive heat of 
formation Ah,°, due to the increased concentration 
of oxygen atoms as 7, rises, more than compensates for 
the negative Ah,,,°. 

The third term is essentially a function of the stagna- 
tion pressure, independent of stagnation enthalpy. The 
enthalpy of the liquid layer, h,, is entirely dependent 
upon the surface temperature, which is in turn a func- 
tion of the stagnation pressure. The magnitude of the 
term [h, — h(s)] is approximately 1,000 B.t.u./Ib. at a 
pressure of 0.1 atm. and 1,600 B.t.u./lb. at a pressure 
of 100.0 atm. 

Thus, to summarize, one can anticipate, and this has 
now been borne out by the correlation of a large number 
of individual calculations, that the effect of increasing 
the flight speed, or equivalently of increasing the stag- 
nation enthalpy is to increase the value of Q*. This 
effect is demonstrated in Fig. 6. The upper curve, 
labeled (Q*)p_19, represents the most optimistic be- 
havior of the material (which occurs during pure sub- 
limation). For a more realistic estimate of the thermal 
response of a refractory oxide, which is subjected to a 
hypersonic environment, Q* must be revised downward 
to allow for the existence of a finite liquid layer. 

Since Q* has been found to be a linear function of IT, 
one may correct the upper curve of Fig. 6 to any other 
arbitrary value of the gasification ratio by writing 


(O*) rer = (O*) par — (1 — T)(AQ*/A4T) = (47) 


where (AQ*/AT), which represents the slope of the 
curve of Q* versus I, has also been correlated, and is 
shown as the lower curve in Fig. G. It is seen that this 
quantity depends only on the flight speed, in much the 
same manner as Q* itself. It is interesting that as 
lr — 0 (pure melting, no vaporization) then 
(Q*)r-0 = (Q*)ra1 — (AQ*/4P) (48) 
which is just the vertical distance between the two 
curves, shown in Fig. 6, at any value of the flight speed. 
In the foregoing analysis we have found it convenient 


to treat [ as an independent parameter. In a certain 
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sense this is realistic, because can in fact be varied 
by altering the viscosity of the liquid phase through the 
addition of contaminants to the solid lattice during the 
fabrication process, which usually act to decrease the 
viscosity. Moreover, we are further justified in adopt- 
ing this technique, since, in so doing, we have evidently 
isolated the effects of stagnation enthalpy from stagna- 
tion pressure. Actually, of course, for a given material, 
I can be determined as a function of the flight condi- 
tions, and, as discussed earlier, the extent of the devia- 
tion of from unity can come either from experiment or 
from theory by matching the gas phase results with a 
liquid phase solution as was done by Scala and Sutton.!” 
Let us therefore proceed with the determination 
of the effects of stagnation pressure. 

In attempting to discern the effects of stagnation 
pressure, examination of Fig. 6 shows that if I’ is con- 
stant, then there does not seem to be any effect of pres- 
sure. However, this does not mean to imply that there 
are no pressure effects. On the contrary, the surface 
temperature 7,, (at which vaporization occurs), rises 
dramatically with stagnation pressure. In fact, it was 
found that the surface temperature of the vaporizing 
ZO: correlates with the stagnation pressure for a given 
value of [. Thus, in the range 0.25 < T < 1.0, the re- 
sults of a large number of calculations indicate that 


Ty, = 4,990 + 530 logiy (P./Po), °R. (49) 


In the range 0 < [I < 0.25, the surface temperature 
undergoes its major variation, and then, at [ = 0, one 
obtains 


T & 4,630 + 460 logy (P./Po), °R. (50) 


The increase in surface temperature with stagnation 
pressure can be explained as follows. At a given value 
of I, Q* is independent of the stagnation pressure, and 
hence the left-hand side of Eq. (16) is independent of 
the pressure. However, at constant I, and a given 
stagnation enthalpy, the heat-transfer rate (the 
numerator) is an increasing function of stagnation pres- 
sure. Since the denominator must also rise in order to 
compensate for the rise in heat transfer—i.e., to main- 
tain a constant ratio—the effective partial pressure of 
the vaporizing species must be larger in order to produce 
the necessary interphase mass transfer. This in turn 
requires a larger surface temperature, which leads ulti- 
mately to the correlation of surface temperature with 
stagnation pressure. Since the viscosity of the molten 
layer drops exponentially with temperature, it is ex- 
pected that the actual I of a vaporizing refractory ma- 
terial which exhibits a liquid phase will decrease as the 
stagnation pressure rises, because the aerodynamic 
shear stresses and pressure gradients will then be quite 
effective in producing convection in the liquid layer, 
and Q* as defined above will decrease. 

Since all of the preceding discussion concerning the 
magnitude of Q* has been based on the assumption of 
€, & O, it is instructive to consider the effects of in- 
cluding the net radiative heat transfer in the analysis. 
It is convenient to introduce Eq. (15) into Eq. (16) 
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by rewriting the left-hand side of the latter as 


on (Q) sno WRp + o VRu (67s! — &vT w') Ja 
= OL) 
Meot. VR, 


Thus, while the left-hand side now depends on the 
parameter éxV Rp, the right-hand side of Eq. (16) is in- 
dependent of this parameter. Also, as noted pre- 
viously, the right-hand side of Eq. (16) is not a strong 
function of the surface temperature. Thus, this form 
of the equation shows that, for a given value of I’, the 
magnitude of Q* must be weakly dependent upon the 
emissivity of the surface. When ézV Rp > 0, the 
numerator and denominator of Eq. (51) decrease 
together. But at constant I’, a decrease in ablation 
rate results in a decrease in interphase mass transfer. 
Because of the exponential dependence of the interphase 
mass transfer upon the surface temperature, the latter 
is reduced slightly as éxW R p > 0, which is entirely 
consistent with the previous statement that Q* is not al- 
tered appreciably. 

The above cause and effect relationship can also be 
explained as follows. Surface reradiation reduces the 
heat which must be blocked or absorbed by the ablat- 
ing material. If the effectiveness of each pound of 
vaporizing material is not altered appreciably, then 
clearly a smaller interphase mass transfer is required, 
which requires a smaller driving surface temperature. 
In general, then, we may conclude that when e¢, VR B> 
0, there is slight decrease in surface temperature ac- 
companied by a sharp decrease in ablation rate. 

Suppose we consider the flight condition shown in 
Fig. 5. Then, for example, upon taking T = 1 and 
be = VR p = 0.245 ft'/?, we find that the surface tem- 
perature is reduced by S80°R., which halves the total 
mass loss. It therefore appears that, in this case, re- 
radiation results in a 2 per cent reduction in surface 
temperature which produces a 50 per cent reduction in 
total mass loss. Upon treating écWVR zp as an inde- 
pendent parameter in the range 0 < éaeWR a = 0.50 
ft.'/? it is found that, for a given value of I, Q* as de- 
fined above does not increase by more than 20 per cent. 
That is, at constant I’, the magnitude of Q* is relatively 
independent of the surface emissivity, although Mot. is, 
of course, significantly reduced. It is also noted that 
ésVR zn > 0 produces an important secondary effect, 
which is, that I itself is altered, somewhat. 

To summarize, then, in this analysis we have found 
that, within the framework of the assumptions made 
earlier, the effective heat of vaporization Q* correlates 
with flight speed (stagnation enthalpy) and gasification 
ratio—i.e., 

Q* = Q*(V., T) (52) 
and since one may anticipate that the gasification ratio 
correlates as 

r =I(V., P,, && Rs) (53) 
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the final result may be written 
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Q* = O*(V,, P., WV Rz) (54) 


and is thus independent of the size of the body when 
€x VR, is small. 


Conclusions 


An analysis has been presented of the heat and mass 
transfer at the surface of a vaporizing refractory oxide 
which was exposed to the environmental conditions en- 
countered during hypersonic flight. 

It was shown that by satisfying the conservation of 
mass and energy at the interface of the vaporizing ma- 
terial, it is possible to determine the mass transfer 
(ablation) rate, the effective heat of vaporization, and 
the surface temperature as a function of the gasification 
ratio I’. 

For the quartz-like material considered, it was found 
that at moderate values of I’ the total mass loss from the 
surface is sufficiently low to make this type of material 
attractive for use as a thermal shield during hypersonic 
re-entry. 

When IT = 0, the effectiveness of the material in ab- 
sorbing heat is merely its heat of melting—i.e., [A,,, — 
h(s) }. 

As anticipated, it was found that the total rate of mass 
loss decreases as I approaches unity (pure sublima- 
tion). This decrease in total mass loss with increase in 
I may be attributed primarily to the contribution of 
two terms. These are, (AQ/Am), defined as an “‘effec- 
tiveness quotient,’ which represents the blocking 
action due to the thickening of the boundary layer, and 
|h, — h(s)] which is related to the energy required to 
break the intermolecular bonds during the vaporization 
process. 

It was also found that the quality Q*, the effective 
heat of vaporization, which includes all heat absorb- 
ing or heat blocking effects at the stagnation point for 
quasi-steady ablation, is an increasing function of stag- 
nation enthalpy, and is independent of body size when 
the ablation process is diffusion controlled (‘‘equi- 
librium”’ vaporization), and the surface emissivity is 

small (e,VR =. = 0.1): 


Appendix 


The compatibility relationship for the element nitro- 
gen can be derived as follows. The condition that the 
surface is impermeable to nitrogen is 


(my + My,)w = 0 (A-1) 
Since the absolute flux of a species is given by the sum 
of the convective and diffusive fluxes—i.e., 

m, = p(v + V,) (A-2) 


one obtains immediately 


| PN Vy + py, = 
Cy + Cy, e 


(A-3) 


My = (edu = iow = 
; = a 
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Upon assuming the applicability of Fick’s law for 
diffusion, the diffusion flux of species 7 is given by 


(piVi)w _ — [pD,,(OC; oy) , - 


”. —py»(Le/Pr)»(OC;/Oy)~ (A-4) 


Thus, in order to evaluate the diffusion flux, one must 
obtain an expression for the concentration gradient of 
species 7 at the surface. Now, the similarity transforma- 
tion utilized in references 7, 14, 18, and 21 requires that 


0/dy = W2( p/p) (due/dx),(0/On)  (A-5) 


However, the similarity between the energy equation 
and the diffusion equation in a frozen gas boundary 
layer means that there exists a unique relationship be- 
tween the temperature gradient and the concentration 
gradient which can be expressed" 


(=) " | Cin ey Cie \() (A-6) 
On/,  LOw — (Cp,/Cp,,) 1 \On/, 


so that, upon combining Eqs. (A-3), (A-4), (A-5), and 
(A-6), one obtains the final result 


Map “ (=) ” 
V 2 pits (dtte dx); Pr}, 
On,{1 — [((Cx + Cn,e/(Cx + Cy,)w] } * 
: P (A-7) 
[4.0 a3 (Cp, Cow) | 
where the eigenvalue 6, is tabulated in reference 21, 
and it is noted that for most surface temperatures of 
interest Cy,, = 0. When this equation is solved simul- 
taneously with Eq. (5) as is shown in Fig. 4, then one 
obtains the unique relationship between the surface 
concentration of the element nitrogen and the surface 
temperature. Upon satisfying the conservation of mass 
and energy it was subsequently found that the actual 
mole fraction of nitrogen at the wall varies between 
limits of 0.2 and 0.5. This indicates that a considerable 
portion of the element nitrogen is driven away from the 
surface during the vaporization process. 
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Wings With Minimum Drag 


Due to Lift in Supersonic Flow’ 


INGEBORG GINZEL* ann HANS MULTHOPP** 


Lhe Martin Company 


Summary 


It has been shown by R. T. Jones! that, in order to produce 
minimum drag, the given lift must be distributed over the wing 
surface in such a way that the sum of the downwash induced by 
this distribution and the downwash induced in reversed flow is 
constant over the wing surface. This combined downwash can be 
expressed by an integral which contains the load as a function 
of the spanwise and chordwise coordinate. The problem of 
finding the appropriate load distribution is thus reduced to the 
problem of finding the solution of a rather cumbersome integral 
equation. 

The severe spanwise singularity of the kernel function is 
handled most easily, as in corresponding subsonic problems, by an 
approximate integration over interpolation polynomials. The 
chordwise load distribution is represented by a limited series 
development in Legendre polynomials. The singularity of the 
kernel function along the Mach lines through any pivotal point 
can be avoided by a similar Legendre development of the com- 
bined induced downwash which is constant. The integral equa- 
tion is thus converted into a system of linear equations for the 
unknown coefficients of the Legendre functions of the load dis- 
tribution at a limited number of spanwise stations. 

Practical calculations are carried out on an electronic com- 
puter. The solutions yield the optimum load distribution and the 
local incidence (twist, camber, etc.) necessary to realize this 
distribution. For many wing plan forms, considerable gains 
over a plane wing appear possible. 


Symbois 
te. = Multhopp coefficients 
b = span of the wing 
c(y) = local chord of the wing 
Le = tip chord 
ey = root chord 
‘5.7? = first integrals of influence coefficients 
ly,»t? = second integrals of influence coefficients 
ty,** = third integrals of influence coefficients 
by," = sth integrals of influence coefficients 
in combined flow 
"A i = sth integrals of influence coefficients in 


two-dimensional flow 


l = Ap/(p/2)U? local load 


m = number of Multhopp stations 
p = pressure 
~ +? = local pressure difference of lower and 
upper surface 
x = chordwise coordinate 
y = spanwise coordinate 


Received November 20, 1958. 
tT The authors wish to acknowledge the valuable help they had 
from Glenn Slike who programed the whole method for the 


electronic computer. 
* Design Specialist, Baltimore Division. 
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R = aspect ratio 
A,i(y), Axl), As(y) coefficients of the Legendre series for 
the angle of incidence 
M = Mach Number 
S wing area 
l incident velocity 
XxX: (Xyt — Xnt)/Cn 
XxX; = Xo + (¢y/en) 
Y (yp — ¥n)B/Cn 
ie = (vy — yn)B/y 
a(x, y) = local angle of incidence 
8 = VM-1 
A = sweep angle of the leading edge 
PAs) Legendre polynomial of order n 
r = downwash angle in combined flow 
y(y) = dimensionless circulation 
u(y) dimensionless moment about the mid- 
point of the section 
Ay) dimensionless moment of inertia 
p density 
: _ Veo 
n cosh~(X/Y) 
Subscripts 
i = induced 
f = forward motion 
r reversed motion 
l = leading edge 
tr = trailing edge 
r = root chord 
c/2 = mid-chord 
vy, n = pivotal sections 
Superscripts 
rs = designing the sth integral of the influence coefficient re- 


sulting from the rth term of the series for the load 


Introduction 


8 bon PROBLEM of flying economically at supersonic 
cruising speed is one of reducing, to as small a value 
as possible, the sum of friction drag, drag due to thick- 
ness, drag due to lift, and interference drag. This 
paper deals with the problem of reducing the drag 
due to lift of the wing alone to its minimum value for 
a given lift. This drag is represented by an area in- 
tegral over the products of local lift and local down- 
wash, which is itself determined by the lift distribu- 
tion. Therefore, the minimum drag must be obtained 
by distributing the total lift in the most appropriate 
way over the surface of the wing. This means that the 
lift distribution is the unknown function of an integral 


equation. R. T. Jones! has given the necessary and 

















sufficient condition that must be satisfied by such dis- 
tribution—namely, that the downwash in the com- 
bined flow field be constant over the plan form—but he 
has not indicated (except for the elliptical wing) how 
the shape of the wing may be determined to satisfy this 
condition. In the present paper, it will be shown how 
this genera! criterion may be employed to determine 
both the lift distribution for minimum drag and the 
actual shape of the lifting surface. 

The nature of supersonic flow leads to a lifting sur- 
face problem. Suitable functions for the load dis- 
tribution in the chordwise direction are used and the 
sum of the downwash induced by forward and rearward 
flow is made constant at certain pivotal spanwise 





sections. 


Derivation of the Integral Equation 


The downwash angle, w/U = aj, is connected to the 
lift distribution 


L = Ap/(p/2)U* (1) 
by the integral 
a(x, y) = —(1/4m) X 


I(x’, y')(x ~ x’ )dx'dy’ 


:... — y’| re y' PV (x aes et i — B(y ia ‘Ff 
(2) 


As pointed out in the literature, the value of this im- 
proper integral depends on the order of integration; 
Eq. (2) presupposes that the integration over x’ is 
carried out first. 

The drag due to lift, D;, is the product of local load 
times local downwash angle integrated over the wing 


area, S: 


Duiol2ye = — 0/9) ITY oss 


I(x, y)l(x', y’)(x — x’) 
= 7 VE -— 27 =FO-7P 


dx'dy'dxdy (3) 


If the integration over x’, y’ is carried out first, the 
condition x — x’ > Bly — y’| says that the integration 
interval consists of the forward Mach cones of all points 
x, yon the wing. If the integration over x, y is carried 
out first, the condition x — x’ > B ly — y’| says that 
the integration over x, y covers the rearward Mach cones 
of all points x’, y’; and the subsequent integration over 
x’, y’ then covers the wing. 

The total lift, Z, is the local load integrated over the 
wing area, S: 


L = (p/2)U2f Sf I(x, y)dxdy (4) 


According to the method of Lagrange multipliers, the 
induced drag for a given total lift is a minimum if the 
expression D; — AL is a minimum for an arbitrary 
constant XA: 
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(5) 


The condition that the first variation with respect to 
the unknown function /(x, y) vanishes leads to the basic 
equation of the problem: 


-1/4e) {| 
x—x’>Bly—y’ 


; I(x’, y’)(x — x’)dx'dy’ 
(y — y’)? V (x — x’)? — By — vy’)? 


(1/47) ff 
x—x’>Bly—y’ 


I(x, y)(x — x’)dxdy 
(y — y’)? V (x — x’)? — By — y’)? 


D; — »L = min 










































—r=0 (6) 


The first term is the downwash angle in the original 
flow [see Eq. (2)]; the second term is the downwash 
angle in reversed flow, so that Eq. (6) can be written 


a,;ta,=X 


This is the theorem of Jones! that for minimum induced 
drag the downwash angle in combined flow is constant 
over the surface of the wing. 

Eq. (6) is an integral equation for /(x’, y’), which 
has some similarity to the integral equation of the sub- 
sonic lifting surface theory; in particular, the singular- 
ity of the kernel function at y > y’ is the same. A 
complete and direct solution of this integral equation i 
appears impossible except for very special wing plan 
forms, such as the ellipses studied by Jones.' 

A function /(x’, y’) is called an approximate solution 
of the integral equation if it fulfills the condition of 
Eq. (6) at a limited number of regularly distributed, 
so-called pivotal stations; and the assumption is made 
that the approximate solution will converge towards the 
real solution with an increasing number of pivotal 
stations. 

The similarity of the integral equation with the sub- 
sonic one would suggest a similar treatment if it were 
possible to overcome the singularity at the Mach line. 
The integrand of Eq. (6) becomes infinite at the Mach 


line if one approaches the Mach line from inside the ir 
Mach cone. Outside the Mach cone, the integrand is Si 
zero. The integral itself remains finite, and it is, there- Se 
fore, possible to overcome the singularity by integra- tl 
tion. tc 

The combined downwash angle is, therefore, repre- A 
sented by a Legendre series: pr 


ay + i > i A,P,(x), P,(x) = l, 


n=1 co 
P(x) = (W12/c)[x — x, — (c/2)], 
P3(x) = (W180/c2){ [x — x1 — (c/2) 2? — (2/12)! (7) 


The approximation consists of keeping the first three 
terms in this series, for which Eq. (6) then yields the 
following condition: 
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MINIMUM DRAG 


- Xtr 
5 1 Xtr I(x’, x — x’)dx'dy’ 
(5) A; = ‘f (ay + a,)Pi(x)dx = — - f ff yk cna dx = 
ct to c ZI 4rc xl (y << y’)? V(x —x $i (y am y bi 
basic Es 1 xe 
A; = - { (ay + a,)P2(x)dx = 
CJn 
ae ff \(x — x’)dx'dy’ ( c 
x— xX, — —}dx = 0 8 
4nc?* owe eS 2 | (8) 
- | Xtr 
A3 = - f (ay + a,) P3(x)dx = 
CJa 
i i ff v’, y’) (x — x’)dx'dy’ ( 4) Pe 
x—-Xx,- _ dx = 
4nci 5 yea (y — y’)? 2 12 
In contrast to Eq. (6), Eq. (8) is written more compactly 
(6) by making the range of integration of the inner double c(yi(x, y) = v(y) + u(y) Po(x) + O(y)P3(x) (9) 
integral cover both the forward and rearward Mach ; ar : : 
: : We can carry as many terms in Eq. (8) as we have in 
iginal cone of point x, y. . Ate ‘ hee 
" the load development, Eq. (9). Eq. (9) is inserted 
was . > The i ies r . 
i P into Eq. (8). The integrations over x’ and x are cz dd 
ritten Definition of the Influence Coefficients 1 (8) ‘sag ante cer 
out, and the resulting integrals are called influence 
The next problem is the nature of the load distribu- coefficients. They are defined by the following nine 
tion itself as a function of x. In the case of the ellipse, equations: 
juced Jones! found a constant load distribution. It is as- i*(y, y’) = 
stant sumed that the best distribution will have no singulari- as Pe ee 
ties at the leading or trailing edges. A reasonably if (x — «') P(x) P,(x)dx'dx (10) 
which general chordwise load distribution is a polynomial. c(y)c(y’) « Vix - x’)? — By — y’)? 
> sub- If this polynomial is orthogonalized in the same way 123 123 
5 2 ‘ "i r,r=i1,24,0 $s = 1, Z, « 
rular- as the series for the combined angle of attack, useful , 
e. A symmetry relations can be expected. Thus, the load The integration interval of the inner integral is again 
ation is written :f the forward and rearward Mach cone of point x, y. 
plan Eq. (8) can then be written: 
lution ai l f= yyy’) + (y, v)uly’) + Fy, vO’) |, 
= — y 
on of 4m y-7F 7 
uted, -19 , ( , -99/ , , -39 , 1) 
seit —— l Ff (y, vy yy’) + 7(y, y’)u(y’) + 7 (y, v)O(y ) ay’ | a1) 
Js the 4dr elt ; 
ivotal . l [a yyy’) + 78 (y, vuln’) + 78(y, »’)O(y’) ly! 
-- dy 
4s yy)? 
> sub- 
were ‘ , : 
; The Influence Coefficients for y ~ y 
b ae , : ios (a) both supersonic and subsonic leading and trailing 
Mach It is advisable to separate the influence coefficients 
7 : é i edges; 
le the into their forward and their rearward values. Con- , ' - 
: 3 oe % ; ee i (b) subsonic leading edge, supersonic trailing edge 
ind is sider the “‘influenced”’ section y and the “influencing iievitedl 
: . a only; < 
there- section y’ of a tapered sweptback wing. Fig. 1 shows : : . ' sien 
ns ‘ ie : iets (c) subsonic leading edge, subsonic trailing edge only. 
tegra- the main possible positions of section y’ with regard , ‘ . ° 
2 Des If y < y’, the situation (a) only occurs. 
to the forward Mach cones of both ends of section y. Br, : ; ; 
- : GE The three cases (a), (b), and (c) are connected with 
repre- At the left of y(y > y’) a greater variety of positions is ote at a i ; 
. . the following integration intervals in the variables x 
possible. 
: bs ; , and x: 
(a) y’ partly inside cone 2 but outside cone 1; 
(b) y’ partly inside cone 2 and partly inside both (a) x) Sx’ Sx — Bly — vy’) 
cones; and x’ +By—y) SxS Xe 
‘ — ols 
c) y’ partly inside cone 1 and completely inside ; = 
3 ( )9 I . I ; (b) "<S x’ Sx -— By —-y’) 
cone 2. a, eae a 
! (7) The situations (a) to (c) are possible for y > y’ and et eee ee (12) 
: (c) x <x’ Sx -— Bly—y’) 
three for x, < x < Xy' + Bly — y’) 
, wv a, 
dis the + The reader is cautioned that the definition of y in reference 2 SxS Xe 


differs from the present one by a factor of 1/4. 








for i" + Bly = y’) < x < Xtr J 
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A similar variety of positions is possible for section cients for rearward Mach cones can be found from 





those for the forward Mach cones. 
As an example, the first influence coefficient i/"(y, y’) 
for torward flow is calculated for the most complicated 


y’ with regard to the rearward Mach cones of both ends 
of section y. It is not necessary to discuss these be- 
cause, as will be seen presently, the influence coeffi- 

case (c) and y > y’: 


j N(y yy") 7 l ztr’ +B(y—y’) x—B(y—y’) (x bean x")dx'dx 
Ay, 9’) = see + 
c(y)c(y’) Je, x’ V(x — x’)? — By — y’)? 
{" : (x — x’)dx'dx l E —%; 
- - - = ~ - Vv (e, — - "2 _— B(y — y’)? 7 
re - ot OCT Ae, ‘ Nir v7 h: ‘ 
Ja +ey—y)J a’ V(x — x’)? — By — y’)? c(y)c(y) 2 
B(y ne y’)? Xe — x,’ x; Cad : , ——— — B?(y — y’)? ; x, x,’ 
- cosh ——_ V (x, — x)? — By — vy’)? + cosh! — 
2 Bly ay oe 2 P “ viet yA Bly — y’) 
Sa Bey — - — By — y’)? , Se Xe e 
s V Cy — %% ) = Po =a) + cosh n (15) 
— a os 2 Bly —- y ) 


As discussed earlier, the influence coefficients are now evaluated at a limited set of m discrete stations. 


y, = sin vr/(m + 1) 


where m is the total number of stations. 


With further abbreviations, 





Y = (y, — Yn)B/Cn 

Xi = [(xX.1 — Xnr)/Cn) + (c,/€n) f= VX2— ¥? m = cosh-(X,/Y) & = = Ofor Y2>X2 

X: = [(x,. -— Xn1)/Cn] ‘ VXx,? — ¥? wm» = cosh-(X./Y) & = m = Ofor ¥2 > Xz? (14) 
X3 = [(X,1 — Xnd/€n] + (G,/en) — 1 & VX;? — Y? 3 = cosh-"(X3/Y) & = n3 = O for Y? > X;? 

the influence coefficients describing the influence of 

station y" = y, on station y = y, for the forward flow ing? = 24/3(c,2/c,2){ (Xi, — Xo%s)/3] — 

can be written V(b, — £:)/3] — (G)/¢a) [Xr — (€,/2cu) Vinny 

bop = (ey Ci (Xie — Xo — X3k3)/2] — Ing? = 6V/5(c,°/c,*) 1 (Xi — X2%fe)/4] — (15) 


(Y?/2)(m — no — 3)¢ (13) 

The expression [Eq. (13)], which holds for the case 
of two subsonic edges, simplifies for other conditions. 
For instance, if the trailing edge is supersonic, the terms 
with index 3 vanish by virtue of Eq. (14); if the 
leading edge is supersonic as well, the terms with 
index 2 vanish by virtue of Eq. (14). For y < y’ 
only the terms with index 1 occur. 


Calculations similar to those in Eq. (13) lead to the 
other influence coefficients.| As an example, we show 
ty)” and 7,,!° for the case of a subsonic leading edge 


and a supersonic trailing edge (& = y; = 0): 


, 


if (y, y’) = 


i 
zU'+B(y 


The integration area 


i [x’ —x/ - (c’ 2) (x — x’ 
-y') Jz V(x — x’)? — B(y — vy’)? 


{" —B(y—y’) , 4 : ; xX) 
J 2)’ J2'+p(y—y’) V (x — x’)? — By — y’)? 


r 


is the shaded area in both cases in Fig. 2. 


These equations hold for » > n. For v < nu (case 
(b), see Fig. 1) only the term with index 1 is used. 

Because of symmetry, it is only necessary to calculate 
for positive values of vy and for the whole range of n. 
Furthermore, one advantage of using the Multhopp 
station is that, for even values of v, only the odd values 
of m are needed, for odd » only the even n. 

Finally, as a result of the orthogonalization of the 
load distribution and of the combined downwash dis- 
tribution, it is not necessary to calculate the rearward 
values of the influence coefficients. We have, for ex- 
ample, for the case of two supersonic edges: 


) 
dx'dx = 


(vw — x’)[x’ — x)’ — (c/2)] 


; —= dxdx'’ = 1,"*(y’,y) (16) 


Similarly : 


T The interested reader can obtain the full set of influence coefficients for any type of plan form by writing to the authors 
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We may, therefore, drop the index f and simply say 


i. = te + ony (18) 


/ 


The Influence Coefficients for y = y 


qs. (13) and (15) are not suitable for the definition 
of the influence coefficients at y = y’ because terms 
like these containing cosh! in Eqs. (13) and (15) lead 
to such as (y — y’)? In (y — y’), andsoon. Inspection 
of the integrals in Eq. (11) shows that the contribution 
of these singularities at y = y’ would remain finite. 
The difficulty is associated with the manner of approxi- 
mating the integrals in Eq. (11) by a finite sum. In 
trying to obviate the difficulty, we are guided by the 
observation that, with increasing Mach Number, the 
integration region in the integrals, such as Eq. (13), 
becomes more and more restricted to the immediate 
neighborhood of station y itself. This suggests that 
one may approximate the local contribution by the 
local contribution of a corresponding two-dimensional 
wing with chord identical to that existing at y. We 
add and subtract the contribution of this fictitious 
wing in such a manner that, on one hand, the dis- 
turbing terms cancel and, on the other hand, they are 
represented by known two-dimensional values. 

For the unswept, untapered infinite wing of chord 
c,, the influence coefficients of Eq. (15) become 


 _—_— ¢* if Y*2y* 
i,'?* = 0 (19) 


_— os V5} Y*?((5/2)&* — n*] — (3/2) Y*'*} 


lyn 
and so on, with 


Y* = (y, — ¥n)B/c,, 


f= Vj — Y*?, »* = cosh—(1/Y*) 
The integrals in Eq. (11) can all be carried out? and 
give in first approximation (as they should) : 
\ = (3/c,) 7, w, = 0, = 0 (20) 
Using now the numerical integration formula’ for 
the integrals of Eq. (11) and adding the two-dimensional 
result of Eq. (20), and subtracting it again in the nu- 
merical integration form, we can write 
2a, = Lo he + Le + "a, —- 
Dean (iru Yn > hag the + hn) 
(Ga,, /c,) aos 4, + hn, + .""a) = 
Lyn (in! * 7, + a + bn "@,) (21) 


The assumption that the local contribution of the 
actual wing can be approximated by the local con- 


t Special care is required in the neighborhood of the wing tips. 
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Fic. 1. The position of section y’ with respect to section y 


tribution of the fictitious two-dimensional wing then 


says 
7 11 — 7 11% 
vy ve 
, 21 — 7 21% 2s 
ly Ps (22) 
5 81 > 31* 
Ly, ly 


which cancels the first three terms on the right-hand 
side of Eq. (21) against the eighth, ninth, and tenth term. 

If we like. we can now write Eq. (21) in the familiar 
form 
2a,,r — i, * "4, + Pigs: rd 

LG en(irn Yn + Fyn™ ttn + Fyn !On) (23) 
by defining new i,,: 
i,,"! = (a,,8/¢,) + layniyn'!* 

i," = 0 (24) 
ion ~ D4 yniyn?!* 


Similarly, we get for the second and third equations 
of Eq. (11) 











Fic. 2. Integration area. 
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Fie = Ln (Tn 2¥n 7 Ta ie + $.—"Ow) 25a) 


ts" Ve + i. "O, = Lyn (Fon? ¥n — ton tn = in *O,) 
(25b) 
with i or (a,,8/c,) + LA nben®?* | 


eo (24) 


= ty 


ie = (a,,8/c,) + LArnin®* | 


The Final Equations 


Eq. (23) mainly determines y (with some interfer- 
ence from yu and 9); Eq. (25a) mainly determines u (with 
interference from y and 9); and Eq. (25b) mainly de- 
termines 90 (with interference from y and yu). As it 
should be, the influence of station ” on station vy de- 
creases with increasing distance between the stations, 
so that the highest terms, for instance, on the right side 
of Eq. (25a), are the neighboring terms 


99 


Byp=itens ’ Getitoss 


99 


and, on the right side of Eq. (25b), the neighboring 


terms 
m4 33 z 33 
Ayy—12tyy—1 » By +itw+tl 


It now turns out that 


are negative, and 7,,”" and 7,,** as defined by Eq. (24) 
are positive. A consequence of this from Eq. (25) is 
that the solutions, u, and 9,, alternate in sign for suc- 
cessive values of v. Such oscillations in u and 0 would 
imply an oscillation of twist and camber which would 
depend on the number of stations. Since this diffi- 
culty is introduced artificially,t the solutions were 
smoothed by an averaging process in which the differ- 
ence between successive » and 0 values is assumed to 
be very much smaller than the absolute local value: 

) 





| — 4, |Ou4, — 9, | 
| My+1 H, | < 1, | : +1 a | < 1 | 
My+1 | 0,41 { ; 
: | ) , (26) 
1- & e,.. — 8, 
| et P<, | tle 
= 0,1 








When one subtracts the two neighboring terms 
(Qy,—12yp—17°* + 4,,417,4172*)u, from both sides of Eq. 
25a) and the two neighboring terms (@,,_17,,—1°°* + 
Ay+thy +1°*)O, from both sides of Eq. (25b), and uses 
Eq. (26), one is led to the final equations: 


io - Dean (inn + Sis ie + 1,n°*O,) (27a) 


ih + t,,**O, = 
Lon! ¥n + Inte + in%O,)  (27b) 


with the definitions 


7 The use of an infinite swept wing rather than of an infinite 
straight wing in deriving Eq. (21) avoids this difficulty but the 
corresponding sum in Eq. (24) then fails to converge. 


tyne? = B97? — 7,,7°* for ly n| = 1} 
in?? = 3,,_%7 for ly — n | > 1 | = 
Zn?? _ Tyn?® a in? 3* for ly _ n | — (28) 
tpn? = 7,4°* for ly — n | >. I | 


and 


i,2? = (a,,B/c,) + do "Qyndon?2* 
Z,,°8 _ (a,,8/c,) 4 > "Gyni,,3** 


where the double omission sign in the sum means that 
the neighboring terms of station v are omitted: 

Eqs. (23), (27a), and (27b) are the final equations. 
They give the constant load distribution found by 
Jones for the ellipses; they give the two-dimensional 
load of Eq. (20) for very high Mach Numbers. 

On the other hand, for an untapered unswept wing 
we get a wing without twist or camber. Graham,® 
using a second-order polynomial in x and y for the 
angle of attack, found some camber for the rectangular 
wing. It is not clear at this stage whether this disagree- 
ment stems from the assumptions of the present theory 
or from the assumptions of Graham, or both. 

The system of Eqs. (23) and (27) divides into two 
systems, one for the odd v and even n, and one for the 
even v and odd n. By inserting one of these systems 
into the other, the number of unknowns can be im- 
mediately divided by two. Using 16 stations over 
half the span, the number of unknowns can be reduced 
to 24. Asin reference 2, the center of a tapered swept- 
back wing must be rounded off slightly. 


The Shape of the Profiles and the Drag 


If the load distribution is known, the shape of the 
profile at section y can be found in terms of the influence 
coefficients i,,; [Eqs. (13) and (15)] of the forward 
Mach cone alone. The terms /,,/” are merely half of 
those of the combined flow [see Eq. (24)]. However, 
superposition of the forward and rearward flows made 
some of the coefficients 7,,’° zero in combined flow be- 
cause of antisymmetry. These now take on finite 
values for forward flow alone. 

The Legendre expansion of the local incidence follows 
the pattern of that of the combined angle [see Eq. (7)], 


a(x, y) = Arly) + Aa(y) Pale) + Aaly)Palx) (29) 
where A,(y) = [1/c(y)] f r a(x, y)dx (30) 
is the twist, 

Axy) = [1/c(y)] f : a(x, 9) Palx)dx 


is called the camber, and 


Aly) = [1/e(y)] f 


Zl 


Ltr 


a(x, y)P3(x)dx 


is called the second-order camber of the wing. 
In terms of the influence coefficients 7,,, in forward 
flow alone twist camber and second-order camber can 
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now be expressed as 


2a,,A 1(y) = is. ae ts = is i. — 
YAyn(in Yn + tyne itn + tyn®!On) 


2a,,A (vy) as ta % aie he + ty es { (a) 
YOu (iyn! Yn = i Sa te + bg Oa) 


20, A 3{ y) aca mt, Y; + e's i i r *¢ My nl 
PAu ( bom Sy, Ss i, othe si i, n ™ Ds ) ) 


which corresponds directly to Eqs. (23) and (25) for 
combined flow and is based on Eq. (2) for the downwash 
angle in forward flow. Aj,(y), Ao(y), A3(y) determine 


minimum induced drag. 
The total lift is 


CL, = & { ; (y/4)dy (32) 

The drag is half the drag in combined flow 
Coi = (A/2) ,. (y/4)dy (33) 
Cp /Cr? = [\/2RS (y/4)dy] (34) 


Discussion of the Results 


The final results depend on three parameters which 
are here chosen to be 


tan A,/B; ci/¢,, ¢,/(0/2) tan A,, 


The first is the ratio of the mean sweep angle of the 
wing to the sweep angle of the Mach line tan A,,/8, 
the second is the taper ratio c,/c,, and the third is the 
ratio of the dimensionless root chord to the mean sweep 
angle of the wing c,/(b/2) tan A,. The mean sweep 
angle was chosen instead of the sweep angle of the 
leading edge because both edges are equally important 
in combined flow. 

Fig. 3 exhibits C,,,,/Cx” for M4 = 3 with aspect 
ratio as abscissa and tan A,, as ordinate. It holds for 
other Mach Numbers as well if the scale of the tan 
Am axis is altered correspondingly. This plot is valid 
for wings with zero taper ratio—namely, for arrow wings 
in the region above the dash-dot curve and for diamond 
wings below this curve. Another curve on the plot 
(dash-dot-dot) separates the subsonic leading edges 
from the supersonic leading edges. Solid curves 
indicate constant values of C,,,,/C,?. The dashed 
curves are those for constant 1/C,, which coincides 
with Cp/C,? for the planar wing if the suction force is 
neglected. The dashed curves, therefore, represent 
upper bounds for Cp/C,? of the planar wing. 

The figure shows that C,,,,/C,? is much smaller 
than 1/C;_, especially for subsonic leading edges. In 
the supersonic leading edge range, the 1/C,, curves 
approach the Cpmin/C,? curves. The mean sweep of 
the diamonds is low; they are, therefore, inferior to the 
arrows of the same aspect ratio as far as minimum drag 
is concerned. 
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Fig. 4 shows the constant Cpmin/C,” curves for wings 
with a finite tip. The curves are steeper than those 
for the pointed wings but they are of the same character. 

Figs. 5 and 6 show pressure distributions of two 
wings at Z = 3, both wings having an aspect ratio of 
1.6 and a sweep angle at the leading edge of 74°. One 
wing is an arrow (Cpmin/Cz,? = 0.5), and the other has a 
finite tip c:/c,; = 1/3 (Cpmin/Ci? = 0.42). Comparison 
of the two figures shows that the character of the pres- 
sure distribution is different towards the tip. For the 
finite tip, the linear term yu of the lift distribution—see 
Eq. (9)—becomes positive and the quadratic term 0 
negative towards the tip (see Fig. 7, which refers to this 
wing). In both cases, the center of pressure is in 
front of the midchord at the root and moves back as 
one proceeds to the tip. 

Fig. 7, which is representative for all calculated cases, 
shows that the third Legendre coefficient 0 of the lift 
distribution becomes small; it can, therefore, be as- 
sumed that the higher coefficients would be even smaller, 
which is a justification for the representation of the 
lift as a Legendre development over chord. Further- 
more, the figure shows that the lift distribution over 
span, y, is almost elliptic. If the optimum load over 
such highly sweptback wings is integrated in spanwise 
direction first and then plotted over the chordwise 
direction, the curve deviates markedly from an ellipse. 
This is in contrast to the load distribution over chord 
of the elliptic wings considered by Jones. 

Since dihedral has no effect on the lift in forward 
flow, an infinite variety of shapes is compatible with 
the theory. The procedure, in fact, gives only the in- 
clination, and an additive constant remains free in the 
equation for the profile at each station. Fig. 8 shows 
two of the possible shapes for the arrow wing of Fig. 5. 
In the upper plot, the wing is twisted about the trailing 
edge; in the lower plot, about the midchord. Fig. 9 
shows in the same manner two possible shapes for the 
wing of Fig. 6. Toward the tip the twist becomes neg- 
ative. The camber is positive everywhere except at 
the root. These shapes were computed for C, = 0.2 
but are also valid for other C, if the scale of the z-axis 
is altered correspondingly. The upturned nose of 
the root section is a characteristic feature of the opti- 
mum wing. 

Generally, one can conclude from the calculated 
cases that the second-order camber, A3;, is small com- 
pared to the first-order camber, A», and it therefore 
can be assumed that the higher order camber terms 
will be small too. 


Conclusions 


The method derives the pressure distribution and the 
distribution of angle of attack for a tapered sweptback 
wing in supersonic flow having minimum drag due to 
lift. Both these distributions are given at 31 stations 
over span in the form of a three-term Legendre series 
over chord. A considerable reduction in drag due to 


(Continued on page 36) 
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Minimum Weight Analysis of Orthotropic 


Plates 


Under Compressive Loading 


GEORGE GERARD* 
New York Unwersity 


Summary 


Minimum weight analysis of various forms of stiffened-plate 


constructior 


multicell wing box beams are presented. 


1 considered for use as the compression covers of 
Longitudinal, trans- 


verse, and waffle-grid stiffening systems for flat plates were in- 
vestigated by use of orthotropic plate theory and the assumption 


that the buckling modes occur simultaneously. 


After obtaining 


results for each stiflening configuration, the comparative struc- 


tural efficiencies of the various types of plates were then con- 


sidered in terms of common parameters. 
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Subscripts 


Ss => 
w = 


x,y = 
1,2,3 = 


Symbols 


exponent 

longitudinal stiffener spacing, in., also exponent 
exponent 

flexural or torsional rigidity, in.-lb. 
modulus of elasticity, psi 

effective modulus, psi 

shear modulus, psi 

moment of inertia, in.* 

torsional constant, in.4 

buckling coefficient 

length, in. 

transverse stiffener spacing, in. 

exponent 

loading, Ib./in. 

ratio of 6 dimensions of stiffener and skin 
ratio of ¢ dimensions of stiffener and skin 
ratio of stiffening system weight to skin weight 
thickness, in. 

effective thickness of stiffened plate, in. 
effective load-carrying thickness, in. 
weight, lbs. 

plate width, in. 

coordinates 

efficiency coefficient of long-stiffened plate 
solidity coefficient 

plasticity reduction factor for plates 
Poisson’s ratio 

solidity 

applied stress, psi 

buckling stress, psi 

buckling stress of orthotropic plate, psi 
local buckling stress of skin, psi 
compressive yield strength, psi 

average or effective stress, psi 

plasticity reduction factor for wide columns 
density, lbs./in.’ 


optimum 

skin 

web or stiffener 
coordinates 
coordinates for D 
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(1) Introduction 


4 tye EVOLUTION of airframe wing structures has 

witnessed a transition from wide-column construc- 
tion, represented principally by stiffened panels, to 
plate construction represented by the multicell or multi- 
web types. This change has been necessitated by the 
thin wings and high loadings associated with the super- 
sonic speed regime. 

A common form of multicell construction utilizes a 
relatively thick cover supported by a series of longi- 
tudinal webs. The basic efficiency of such a structural 
arrangement in thin wings has encouraged investiga- 
tions concerned with the rearrangement and reduction 
of structural material within the cover plate with the 
objective of obtaining decreased solidity. 

It is important to note that the covers of the wing 
perform a contouring as well as a load-carrying func- 
tion. This fact together with stability or buckling con- 
siderations can result in an optimum stress for the com- 
pression cover which may be considerably below the 
compressive yield strength for specified loadings and 
supporting structure configurations. 

The approach used herein is to analyze idealized 
long-plate structures which retain the significant de- 
tails of the actual structure and yet are sufficiently sim- 
plified so as to permit broad conclusions to be drawn 
as to the optimum stress and configuration of the 
minimum weight plate. A comparison of the optimum 
plate of each type in terms of common weight and load- 
ing parameters can then be etfected to establish the 
ranges of efficient application of each type. This is es- 
sentially the approach required in design, as distinct 
from stress analysis, in that the minimum weight 
structure for a given set of loading and geometric 
parameters is desired. 

The analytical techniques which are available for the 
minimum weight analysis of structural elements exist 
primarily for wide-column or stiffened-panel types of 
construction. Such analyses have been presented by 
Zahorski,! Farrar,? and Catchpole.* A comparison of 
the results of such analyses with experimental data is 
contained in reference 4, in which good agreement is 
evident. 


(2) Long Orthotropic Plates in Compression 


In considering various types of stiffened plates, it is 
pertinent to denote the type according to the orienta- 
tion of the stiffening system. The major stiffening sys- 
tems for single-skin construction are illustrated in Fig. 1 
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Fic. 1. Classification of stiffening systems for single-skin con- 
struction: (a) longitudinal, (b) transverse, (c) grid (longitudinal 
plus transverse). 


and include longitudinal, transverse, and grid systems. 
As shown in Fig. 1, the grid system is essentially a com- 
bination of longitudinal and transverse stiffeners, al- 
though grids at various angles to the plate axes can be 
included in this category. 

When the stiffener spacings for the plates shown in 
Fig. 1 are small relative to the overall length and width 
of the plate, it is possible to use orthotropic plate 
theory to determine buckling stresses. In this case, the 
plate has different rigidities in the two orthogonal direc- 
tions in its plane, and because of the closeness of the 
stiffening system the actual plate may be represented 
by an equivalent uniform orthotropic plate. 

Timoshenko’ has presented the buckling analysis of a 
long, simply supported, flat orthotropic plate under 
axial compression. The governing equilibrium equation 


1S 


D,(04w/Ox*) + 2D;(04w/Ox?dy?) + 
D.(0*w/dy*) + N,(0?w/dx?) = 0 (1) 
where 
dD, (EI),/(1 — vzvy); average flexural rigidity of 
a unit transverse cross section 
D, = (EI),/(1 — »v,v,); average flexural rigidity of 
a unit longitudinal cross section 

(1/2)(v,D2 + vyD,) + 2(GI)r7,; where (G/),y is 

the average unit torsional rigidity 


Il 


Ds; 


Upon substitution of a suitable solution for Eq. (1), 
the following equation for compressive buckling of a 
long orthotropic plate is obtained: 


Ocr = (22?/wt)[(D,D2)'” + Ds] (2) 


For an isotropic plate, D; = D, = D3, and Eq. (2) re- 
duces to the familiar result for this case. 

The flexural and torsional rigidities of a stiffened 
plate can be calculated if the plate is monolithic and 
the stiffening system is of simple configuration. For 
built-up construction and more complex configurations, 
experimental methods are generally employed. 

The theory based on Eq. (1) applies to stiffening sys- 
tems which are symmetrical with respect to the middle 
surface of the skin. Asymmetry of the stiffening system 
results in the centroid of the section not being in the 
plane of the skin and thus introduces additional terms 
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in Eq. (1) for such cases. Dow, Libove, and Hubka® 
have presented the more exact theory for single-skin 
construction. 

The asymmetry of the stiffening system causes bend- 
ing to proceed simultaneously with axial loading, and 
thus the results of a buckling analysis for this case re- 
semble those for an initially imperfect column. Craw- 
ford’ has presented the results for an integrally stiffened 
plate with a one-sided stiffening system. The results 
indicate that, although some small amount of bending 
is evident immediately upon the application of axial 
load, the maximum load-carrying ability of the plate is 
essentially the same as the critical load obtained from 
Eq. (2). 


(3) Longitudinally Stiffened Plates 


In this section, the minimum weight analysis of a long 
longitudinally stiffened plate which is assumed to be 
simply supported along the unloaded edges is con- 
sidered. The characteristic form of this plate is shown 
in Fig. 1, and it is assumed for the purposes of analysis 
that the cross section is as shown therein. The use of 
an unflanged monolithic stiffening system is of con- 
siderable significance in thin-wing construction, quite 
aside from the simplicity of its form from an analytical 


standpoint. 
Geometric Properties 
The effective thickness of the section shown in Fig. 
1(a) is given by 
i/t, = 1 + Trl t (3) 
The average flexural rigidity of the cross section parallel 
to the loaded edges can be obtained in the following 


form 


Dy = [nE, (i -— Vrvy) \(dwts ‘12) x 
rer 1[(4 + Yo? 1) (1 + Yor 1) | (4) 


Eq. (4) is valid for many closely spaced stiffeners. Ina 
practical sense, it applies for three or more stiffeners. 

The flexural rigidity perpendicular to the loaded 
edges is simply that of the skin alone. 


D, = D, = nF t,3/(12(1 — v,,)] (5) 


The torsional stiffness of the plate, D;, is taken to be 
small since, by virtue of the reciprocity relation, 


Vy ‘Dy = Ds (6) 


and, by use of the relationship, 


“J 


J=4i,, ( 
the definition of D; from Eq. (1) becomes 
D3 = vzDe + (GJ 2) (8) 


Since the torsional rigidity of the stiffeners is small, it 
follows from Eq. (8) that D; must be of the same order 
of magnitude as D2, which itself is small relative to D,. 
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Buckling Modes 

The assumed conditions for minimum weight design 
of a compression structure are that the applied stress 
and allowable stress be equal and that the possible forms 
of buckling occur simultaneously. The applied stress in 
this case is 

a. = N/i (9) 

Overall buckling of the long, simply supported, longi- 
tudinally stiffened plate is assumed to be given by the 
orthotropic plate-buckling formula, Eq. (2), with the 
additional assumption that D; is negligible as compared 
to the other rigidity terms. Thus, Eq. (2) becomes in 
this case, 


(ocr) » = (24?/w*l)(D\D,)” (10) 


The stiffened plate is then considered as an assembly 
of plate elements, and a second type of compressive in- 
stability involving local buckling of these elements is 
given by 

(Ger), = kyw*D,/b,*ts (11) 

The general objective at this point is to combine 
Eqs. (9)—(11) in such a manner that the resulting stress 
equation contains only the geometric terms 7, and r, and 
the loading parameter pertinent to this problem, 
N/w. For the stiffened plate, the procedure is some- 
what more involved than for wide columns (reference 2) 
and, therefore, the following approach is required. 

It can be observed that, exclusive of the 7, and ¢, terms 
and certain constants, from Eqs. (9) and (3) 


ag ™~ N/t; (12) 
From Eqs. (10), (4), and (5) 
(Ger) » ~ Dut ,/w? (13) 
From Eqs. (11) and (5) 
(ocr), ~ (ts/bs)? (14) 
By combining Eqs. (12)—(14) in the following form, 
(Ger) ,"(Ger) p'(oa)° ~ (ts/bs)°*(Buts/w*)’(N/ts)° (15) 
Upon the conditions that the margins of safety are zero, 


(ser)p — oe = 04 (16) 
(ocr); — Oa = of 
By use of dimensional analysis, the exponents take on 
the following values: 


= 1 
nin | (17) 
c= 44 


and, denoting the stress under the conditions of Eqs. 
(16) as @, Eq. (15) becomes 

a’ ~ 1,?(N/w)! (18) 
This relation represents the fulfillment of the desired 


objective. 
By using this same procedure for the complete Eqs. 


(9)—(11) and the exponents given by Eq. (17), the fol- 
lowing stress equation is obtained for the longitudinally 
stiffened plate with unflanged monolithic stiffeners. 


& = a,[fE/1 — v,v,)|/C"(N/w)" (19) 
where 
a, = Fal " [hater (4 + ror) . (20) 
12° 1 + rr, 
i = (2m)? (21) 


Eq. (19) is in a form similar to that presented in ref- 
erences 1-4 for wide columns. It expresses the de- 
pendence of stress upon the loading index, N/w, the 
effective modulus and the stiffened-plate efficiency co- 
efficient, a,. The latter is solely a function of the 
geometric ratios 7, and r, and buckling coefficient k, 
and as such can be maximized to obtain the greatest 
efficiency for a prescribed cross-sectional configuration. 


Optimum Longitudinally Stiffened Plates 


Exact methods of analysis for optimizing the ef- 
ficiency coefficient have been presented in some detail in 
reference 2. These methods were developed for wide 
columns and apply directly to stiffened plates. For the 
purposes of this analysis, it is sufficiently accurate to use 
the approximation that the various elements in the 


cross section are pin-jointed. For simultaneous 
buckling of the plate and stiffeners, then, 
0.43(tw/bu)? = 4(t,/bs)? (22) 
therefore, 
r, = 3.05 7, (23) 


By substituting Eq. (23) into Eq. (20), the stiffened- 
plate efficiency coefficient becomes, for k, = 4, 


ees ‘7 (4r,4 + 3.05r,%)!/7 
a, = 


(24) 
125 1 + 3.05r,? 


Upon maximizing Eq. (24), 0a,/0r, = 0, it is found 
that 


(r,), = 0.376 (25) 
By substituting Eq. (25) into Eq. (20), 

(a;), = 0.796 (26) 
Finally, by substituting Eq. (25) into Eq. (23), 

(r,), = 1.146 (27) 


It can be noted from Eq. (3) that the weight of the 
stiffening system relative to the weight of the skin is 
S = ref; (28) 
For the optimum longitudinally stiffened panel, from 
Eqs. (25) and (27), then, 
S, = 0.431 (29) 
By virtue of Eq. (26), Eq. (19) becomes 


&, = 0.796[7E/(1 — vzvy) }*7(N/w) (30) 
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In comparing the efficiency of various types of con- 
struction, it is best to make the comparison on a weight 
basis rather than a stress basis, since the latter often 
does not include the weight of the supporting structure 
(reference 4). For this purpose, the weight of a planar 
structural element is given by 

W = glwL (31) 


The solidity or structural density of this element can be 
defined as 
> = W/ew*L = i/w = (N/w)/e (32) 
By substituting the optimum stress relation, Eq. (30), 
into Eq. (32), 
>, = 1.257{[(1 — vv) N]/qEw}*” (33) 


It can be noted that Eq. (33) is in nondimensional form. 


(4) Transversely Stiffened Plates 


The minimum weight analysis of a long transversely 
stiffened plate, such as shown in Fig. 1, simply sup- 
ported along the unloaded edges is essentially the same 
as for the longitudinally stiffened plate of the previous 
section. Again, an unflanged monolithic stiffening sys- 
tem is assumed which in this case acts to stabilize the 
skin but does not perform any load-carrying function 
as does the longitudinal stiffening system. In such a 
case, it is necessary to distinguish between an effective 
thickness for load carrying and one from a weight stand- 


point. 


Geometric Properties 
The effective thickness of a longitudinal cross section 
from a weight standpoint is 
i/tj =1+ Hr, (34) 
where in this case, 7, = b,/L,, and L, is the spacing of 
the transverse stiffeners. Parallel to the loaded edges, 
the flexural rigidity is that of the skin alone 
D, = D, = 7,Et,*/(12(1 — very) ] (35) 
Perpendicular to the loaded edges, the average flexural 
rigidity is 
Dz = [nE/(1 — VzVy) |(dw7ts/12) X 
ror, (4 + rors)/(1 + ror.)] (86) 
By following the same argument given in Section (3), the 
torsional rigidity, D3, is assumed to be negligible. 
Buckling Modes 


Since the stiffening system does not carry any axial 
load in this case, the applied stress is 


og = N/t, (37) 

Overall buckling of the transversely stiffened plate is 
given by 

(Ger)y = (24?/w*t,)(D,D2)'” (38) 


The local buckling mode involves the skin as a wide 
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column simply supported at the transverse stiffeners 
since the latter have been assumed to have negligible 
torsional rigidity. Thus, 
(Ger), = wD,/L/*t; (39) 
By following the identical procedure used in Section 
(3), Eqs. (37)—(39) can be combined to obtain the follow- 
ing stress equation: 


& = a,[7E/(1 — v,v,) |*/"(N/w)”” (40) 
where 
Bo ty (i 1/7 
ai E l r(4+ 2] (41) 
123 (14+ nr, 


j= (7,2n)”” (42) 


Since the load-carrying structure is less than the total 
structure in this case, it is necessary to proceed directly 
to the weight or solidity equations before optimizing the 
cross-sectional dimensions. For the transversely stiff- 
ened plate, 


= i/w = [(N/w)/e)(i/ts) (43) 
By substituting Eqs. (34) and (40) into Eq. (48), 
DL = BE — vary) N]/aEw}*” (44) 
where 
B = (1+ nmr,)/a, (45) 


Efficient Transversely Stiffened Plates 

In a transversely stiffened plate, the stiffening sys- 
tem performs only a supporting function and not a load- 
carrying function. Since the stiffener is not under axial 
load, it is not subject to buckling, and therefore there 
is no condition to provide the needed relation between 7, 
and 7, as in the case considered previously. 

To determine the conditions for efficient design of a 
transversely stiffened plate, then, it is possible to mini- 
mize Eq. (44) in terms of the supporting structure 
weight 

S = nr; (46) 


The value of S is determined from Eq. (34). 
By performing the operation, 02/0S = 0 on Eq. 
(44), it is found that 
S, = 0.15 (47) 
Consequently, 
(r:), = 0.15/1% (48) 
By substituting into Eq. (45), 
B, = 1.12/(r,)*" (49) 
It is interesting to observe that Eq. (44) for the 
transversely stiffened plate is identical in form with 
Eq. (33) for the longitudinally stiffened plate. In fact, 
for 7, = 0.67, 8, = 1.257 in Eq. (49) and the efficiencies 
of both stiffening systems are identical, although the 
proportions are different. 
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(5) Grid-Stiffened Plates 


The grid form of stiffening system shown in Fig. 1 is 
essentially a combination of the longitudinal- and trans- 
verse-stiffening systems considered in the previous sec- 
tions. The transverse stiffeners act to reduce the effec- 
tive length of the longitudinals but have little influence 
upon the buckling stress of the skin. This is due to the 
fact that the longitudinals force the skin to act as a 
series of plates which are little influenced by length. 

The minimum weight analysis of a long, grid-stiffened 
plate simply supported along the unloaded edges is pre- 
sented for a particularly simple grid configuration. A 
square-grid arrangement of the type shown in Fig. 1 is 
used in an attempt to assess the comparative efficiency 
of this type of construction. A review of reference 6, for 
example, reveals that a wide variety of grid configura- 
tions is possible, although test data of reference 8 ap- 
parently reduce the number of configurations that are of 
practical importance. 


Geometric Properties 

The grid-stiffened plate to be analyzed has a square 
0°-90° grid configuration with longitudinal and trans- 
verse stiffeners of identical shape and spacing, L, = ),. 
The effective thickness of a unit element of the stiffened 
plate from a weight standpoint is 


t/t, = 1 + 2nr, (50) 
Since the transverse stiffeners do not participate in 


carrying the axial load, the effective load-carrying thick- 


ness is 
A/tp=lt on | (51) 


Because of the symmetry of the grid configuration, 
the flexural rigidities D, and D» are identical. 


Dz. = D, = [nE/(1 — v?)](bn7t,/12) X 

ror [(4 + YP 1) (1 _ Yr? 1) | (52) 
As with the other single-skin types of construction, the 
twisting rigidity D; is assumed to be small as compared 
to dD, 


Buckling Modes 
The applied stress for the grid-stiffened plate is 
oq = N/ta (53) 
Overall buckling is given by 
22°D,/w'i, (54) 


(Ger) » = 


Local buckling involves the skin as a square plate which 
is assumed to be simply supported along four sides due 
to the negligible torsional rigidity of the stiffeners. 


Thus, 
(cer), = {4a E/[12(1 — v?)]} (t,/b,)? (55) 


By following a procedure similar to that used in Sec- 
tion 3, Eqs. (53)—(55) can be combined to obtain the 
following stress equation: 
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& = a,("E/(1 — v?)]'?(N/w)'” (56) 
where 
e ( 8 ) * Iry®r (4 + ror) 4 = 
a,=f (O74) 
beg 1+ nmr, 
i = (nn)'”” (58) 
The solidity of the grid-stiffened plate is 
> = i/w = [(N/w)/e](7/,) (59) 
By substituting Eqs. (50), (51), and (56) into (59), 
> = Bf [1 — »)N)/aEw)'” (60) 
where 
B = (1 + 2nr,)/[as(1 + rr.) ] (61) 


Optimum Grid-Stiffened Plates 


For the grid-stiffened plate under consideration, the 
transverse stiffening system is identical with the longi- 
tudinal system. By assuming that the various ele- 
ments in the cross section are pin-jointed at the line of 
attachment to the skin, it follows that for equal buckling 


stresses for the skin and longitudinal stiffeners 
7, = 3.051, (62) 
By substituting Eq. (62) into Eq. (61), 
9\l/2 6 2 
tear (12)""(1 + 6.177) | (63) 
rr, (24.4(4 + 3.05r,”) }'/* 


By minimizing Eq. (63), 08/07, = 0, the following re- 


sults are obtained: 


(r,), = 0.425 (64) 
(r,), = 1.30 (65) 
B, = 1.675 (66) 
(a;), = 0.81 (67) 
Ss. = 1.106 (68) 


The results given by Eqs. (64)—(68) are for the 0°-90° 
grid-stiffening system with identical and equally spaced 
longitudinal and transverse stiffeners. More generally, 
however, Eqs. (56) and (60) are valid for any grid- 
stiffening system in which D, = D.. The values of a, 
and 8 will, of course, differ from those given by Eqs. 
(64)—(68). 

Test data of reference 8 indicate that of several grid 
configurations investigated, the 0°-90° system such as 
analyzed herein is quite efficient. In fact, only the 
+45° system is somewhat more efficient and appears 
to be the best practical design to carry both compres- 
sion and shear loads. For the +45° system, D,; = Dz, 
and therefore Eqs. (56) and (60) for ¢ and &, 
The +45° system, however, possesses 


respec- 
tively, apply. 
greater efficiency by virtue of the facts that both sets of 
stiffeners share in carrying the load and the buckling 
stress of the skin is increased considerably. 

The test data indicate that the +45° system is ap- 
proximately 50 per cent more efficient than the 0°-90 
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Fic. 2. Comparative structural efficiencies for compressive 


elastic buckling of long, simply supported plates. 


system. Therefore, as an estimate of the optimum 
+45° system, the following value is used: 


B, = 1.12 (69) 


(6) Comparative Structural Efficiencies 


The minimum weight analysis of various types of 
stiffened plates have been considered in the previous 
sections. The specific configurations studied are repre- 
sentative of the types of construction that have been 
proposed for current designs, although the analyses are 
necessarily based on various assumptions and idealiza- 
tions. In this section, the assumptions and results for 
the different types of construction are summarized for 
the purpose of investigating the comparative structural 
efficiencies of the various configurations. 


Summary of Results 

The optimum solidities of the various forms of stiff- 
ened plates all have a common form represented by the 
following: 


d. = B,(N/Ew)” (70) 


The £ term in Eq. (70) represents the plasticity reduc- 
tion factor and Poisson ratio terms appearing in each 
of the equations. In a generalized form, 


E = 4E/(1 — »v,) (71) 


The analytical results for the various types of con- 
struction are summarized in Table 1, together with the 
major assumptions involved. It can be observed that 
the stiffening systems can be classified according to the 
value of m appearing in Eq. (70). This fact immedi- 
ately permits an investigation of the relative structural 
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efficiencies of the forms of construction with the same 


value of m by noting the relative values of 8,. 


Comparative Efficiencies 

For the m = 3/7 case, longitudinally and trans- 
versely stiffened, long, simply supported plates are of 
comparable efficiencies according to the respective 8, 
values. A small increase in efficiency can be obtained 
by increasing the 7, value for the transversely stiffened 
plate although this degree of freedom is not permitted 
for the longitudinally stiffened plate. 

In order to investigate the relative efficiencies of the 
various m categories, it is necessary to resort to a plot 
of the results. Before doing so, however, it is of value 
to consider the buckling of a long, simply supported, 
unstiffened flat plate subject to compression, for refer- 
ence purposes. 

The optimum stress for buckling of a long, simply 
supported, unstiffened plate can be derived from the 
usual buckling equation in the following form: 

2/3 


a, = 1.46[n,£/(1 — v*)}3(N/w) (72) 


In terms of optimum solidity, 
>, = 0.685{ [1 — v?)N]/n,.Ew} 8 (73) 


It can be observed that Eq. (73) is in the general form 
of Eq. (70). 

The comparative structural 
various forms of stiffened plates listed in Table 1, as 
well as unstiffened plates, are shown in Fig. 2. The 
sloping lines toward the left of this figure represent 
elastic buckling, whereas the 45° lines toward the right 
represent strength limitations as opposed to stability. 
The 45° strength lines are drawn for values of E/o,y of 
100 and 200, which encompass a range pertinent to 
many materials at both room and elevated tempera- 


efficiencies of the 


tures. 


Conclusions 

On the basis of this analysis and for the various as- 
sumptions involved, the following conclusions may be 
drawn from the elastic buckling data presented in 
Fig. 2: 

(a) Longitudinally and transversely stiffened plates 
in compression represent only a small improvement in 
efficiency over unstiffened plates. 

(b) Both 0°-90° and +45° waffle-grid construction 
represent a significant increase in efficiency over the 
unstiffened plate. 


(Continued on page 64) 











TABLE 1 
Summary of Results on Optimum Forms of Long, Simply Supported, Stiffened Plates 
Plate Configuration Bo m So Assumptions 
a) Longitudinally stiffened 1.26 3/7 0.43 Unflanged stiffeners of zero torsional rigidity; 
pin jointed 
b) Transversely stiffened ro/0.67 1.26 3/7 0.15 Unflanged stiffeners of zero torsional rigidity 
1.00 1.12 3/7 0.15 
c) Grid-stiffened 0°-90° 1.68 1/2 rot Unflanged stiffeners of zero torsional rigidity 
+45° 1.12 1/2 Estimate based on 0°-90° analysis plus test data 
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Turbulent Heat Transfer on Blunt-Nosed 
Bodies in Two-Dimensional and General 
Three-Dimensional Hypersonic Flow! 


ROBERTO VAGLIO-LAURIN* 
Polytechnic Institute of Brooklyn 





Summary Yo = reference length 
R = mean density 
Recent results obtained for three-dimensional laminar bound- Re = Reynolds Number [Eq. (25) 
ary layers are extended to the turbulent case. It is shown that t = time 
in the presence of highly cooled surfaces and of moderate Mach T = temperature 
Numbers of the outer stream, the crossflow and the pertaining VN), V2, v = fluctuation velocity components 
Reynolds stresses in a general three-dimensional turbulent Vi, V2, Vs; = mean velocity components 
boundary layer are negligible even for large transverse pressure Vie* = (r*/po)”? friction velocity for the transformed 
gradients. A correlation due to Mager between two-dimen- flow 
sional compressible and incompressible turbulent boundary 1, Xo, X = orthogonal curvilinear coordinates 
layers is extended to the problem in question. From a study of 5 = dimensionless variable {[Eq. (21c)] 
the transformation and of its implications, a rapid method for 3 = boundarv-laver thickness 
the analysis of the boundary-layer flow under the subject con- ’ = momentum thickness 
ditions is established. In the absence of general three-dimen- « = component of streamline curvature in plane 
sional data, a comparison with experiments and with the pre- tangent to surface 
dictions of other known analyses is carried out for several axi- ue = viscosity coefficient 
symmetric configurations; the results of the method presented y = kinematic viscosity coefficient 
here exhibit good agreement with the data. The range of valid- ti, ¢ = transformed curvilinear coordinates [Eq. (11) 
ity of the cold wall approximation for general three-dimensional p = fluctuation of density 
problems is estimated qualitatively on the basis of recent meas- “ = shearing stress tensor 
urements in laminar flow, the argument being that, for either Te = shearing stress at the wall 
zero or favorable streamwise pressure gradients, smaller three- y = stream function [Eq. (12) 
dimensional effects are to be expected in a turbulent boundary 
layer, as compared to a laminar layer. Subscripts 
‘ e = at the outer edge of the boundary 
Symbols 0 = reference state used in the transformation 
ct = local skin-friction coefficient r = average oe conditions at a given station 
Cp = coefficient of specific heat at constant pressure oat _—— pnd 
i cae . ge a s = stagnation conditions of the outer flow 
Ci, ¢ = quantities defined by Eq. (21b) ' f 
D = quantity defined by Eq. (21c) ‘i = Slee 
def V = deformation tensor 5 c 
“1 . a B 
ey, C2, € = length elements for curvilinear coordinate mperscripes 
system ‘ fluctuations 
h = enthalpy ° = transformed conditions 
H = stagnation enthalpy indicates time average 
sn = stagnation enthalpy ratio |Eq. (14a)] 
k = thermal conductivity = 
K = proportionality constant in the law of the wall (1) Introduction 
[Eq. (20a)] : ‘ . : 
HE GAP between the available information con- 
K = heat conduction vector ‘ : : 
Vn = Nusselt Number [Eq. (25)] cerning the analysis of laminar boundary layers 
p = fluctuation of pressure and that pertaining to turbulent layers is increasing 
P = mean pressure with the speed of flight. Thus, while problems in two- 
> = Jerr Y, . ° e . *,* 
l Prandt] Number «ol dimensional and, at least under certain conditions, 
Ou = heat transfer per unit area per unit time ‘ . 3 
three-dimensional compressible laminar boundary lay- 
5 ; ‘ ers can be considered well in hand, the recent advances 
Received December 8, 1958. Revised and received May 10, : F 
1959 in the understanding of turbulent flows, as exemplified 
Joe. 


by the work of Clauser,! Coles,? and Townsend,’ are 
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still connected with the two-dimensional incompressible 


tract AF 33(616)-3265 monitored by the Aeronautical Re- 


search Laboratory, WADC, USAF. The author wishes to ex- case. Even under the latter conditions, the methods— 
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face the problems associated with the hypersonic flight 
of missiles and aircraft. The determination of the heat 
transfer on bodies and wings acquires particular in- 
terest in the hypersonic speed range, with laminar, 
transitional, and fully developed turbulent flows being 
all of practical concern. The large rates of heat 
transfer and the strong pressure gradients to be en- 
countered under the aforementioned conditions are 
quite remote from those considered by the semiem- 
pirical incompressible turbulent boundary layer theo- 
ries; hence, simple extensions of classical low-speed 
methods to the new range of interest are debatable even 
for two-dimensional flows. This point of view is evi- 
denced by a critical evaluation of such analyses per- 
formed by Libby and Cresci.!* 

Correlations between two-dimensional compressible 
and incompressible turbulent boundary layers on in- 
sulated surfaces have recently been investigated by sev- 
eral authors; satisfactory agreement with experiments 
is obtained thereby. The analyses of Englert’ and 
Mager® based on the integral and on the differential 
equations of motion, respectively, are typical of this 
group. The extension of Mager’s approach to the 
hypersonic case (in the presence of pressure gradient 
and heat transfer at the wall) is included in the present 
paper 

The problem of three-dimensional turbulent bound- 
ary layers has received little attention even for low- 
speed flows, mainly because the pertaining turbulent 
stress representation is not well established. Very 
little information is available. Mager’ developed an 
integral method whereby the growth of the boundary 
layer is analyzed following the streamlines of the in- 
viscid flow and assuming that skin friction and other 
properties satisfy the known semiempirical two-dimen- 
sional laws. Braun® has recently extended this ap- 
proach to compressible flows over arbitrary curved sur- 
faces. Due to the uncertainty in the selection of proper 
velocity profiles for the crossflow, the quantitative 
accuracy of the integral method predictions is open to 
question. However, the qualitative conclusion ob- 
tained is that, in the presence of favorable streamwise 
pressure gradients, significantly smaller crossflows and 
concomitant effects are to be expected in a turbulent 
boundary layer, as compared to a laminar layer. On 
physical grounds this behavior can be attributed to the 
higher shear stresses encountered in the turbulent case. 

Of particular current interest is the determination of 
the heat transfer on axisymmetric blunt-nosed bodies 
at incidence and on sweptback wings with blunt edges 
in hypersonic flow. The laminar boundary layer on the 
nose portion of general three-dimensional blunt bodies 
in hypersonic flight has been investigated in reference 
9, where it has been shown that the small crossflow 
approximation is valid for an arbitrary streamline 
pattern in the presence of highly cooled surfaces and of 
moderate Mach Number outer flows; hence, two- 
dimensional methods of analysis can readily be extended 
to the three-dimensional problem. The validity of the 
two-dimensional approximation for laminar flow is 


justified by observing that the momentum equation in 
the direction normal to the outer streamline is homo- 
geneous when the relation (R,/R) = (Vi/ V1.) is satis- 
fied across the boundary layer, X being the density and 
V, the velocity component in the streamline direction. 
The physical justification lies in the fact that when 
Ry/R, > 1, the velocity and, even more, the enthalpy 
profile near the surface are affected only slightly by 
the pressure gradient, in contrast to the more familiar 
case of moderate temperature differences across the 
boundary layer. 

In the present paper, an analogous conclusion is ob- 
tained for turbulent boundary layers under the subject 
hypersonic conditions by establishing two equations 
homogeneous in the crossflow and in the pertaining 
Reynolds stresses. Again the demonstration is based 
on the relation (,/R) = (Vi/V;,.)?F and thus requires 
that Reynolds analogy be satisfied by the flow in the 
The validity of this approxi- 
extending 


streamline direction. 
mation is established 
Mager’s correlation® between compressible and incom- 
pressible turbulent boundary layers. From a study of 
the transformation a rapid method for the analysis of 
the boundary layer under the subject conditions is 
developed. Due to unavailability of accurate pres- 
sure and turbulent heat-transfer data for blunt bodies 
in three-dimensional flow, a comparison is instituted 
with the measurements made by Libby and Cresci!*® '4 
on three axisymmetric configurations. Fair agree- 
ment with the experimental data and an improvement 
on the predictions of other known analyses are obtained 
by the method of this paper. 

The various portions of the investigation are de- 
scribed in the following sections according to the se- 
in the previous paragraph. The 


subsequently by 


quence outlined 
utilization of pressure measurements and of theoretical 
analyses in determining the inviscid streamline pattern 
on a three-dimensional blunt body is discussed in refer- 
ences 9 and 10. 


(2) Analytical Considerations for Three- 
Dimensional Boundary Layers 


With reference to an orthogonal curvilinear coordi- 
nate system (x1, x2, x3) having the surface .; = 0 coin- 
cident with the surface of the body, the equations 
governing the steady mean motion in the turbulent 
boundary-layer flow of a homogeneous gas, in the ab- 
sence of body forces and of heat sources, are 


(0/Ox}) [evesRVi] + (0/Ox2) [ere3R Vo] + 
(0/Ox3) [e1e20(RV3 + pv3)] = O (la) 
Vi(OV1/e:0x1) + V2(O0V1/e20x2) + [V3 + (pv3/R)] X 
(OV;/e30x3) + (Vi V2/exe2)(0e1/Ox2) — 
( V22/e1e2) (Oe2/0xX1) = —(1/R)(OP/e,0x;) + 
(1/R)(0/e30x3) [u(OV1/e3:0x3) — Rzyw3] (1b) 


1 R and V; represent mean flow properties in the turbulent 


case, 
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V(OV2/e0%1) + V2(OV2/e20x2) + [V3 + (pv3/R)] X 
(OV2/e30x3) + (Vi V2/e1e2)(Oe2/0x1) — 
(V1?/eye2) (Oe,;/Ox2) = —(1/R)(OP/e20x2) + 
(1/R)(0/e30x3) [u(OV2/e30x3) — Rvevs3] (1c) 


P = P, + Ofs] (1d) 


Vi(OH/e,0x1) + V2(0H/e.0x:) + [Vs + (pv5/R)] X 
(OH /e30x3) = (1/R)(0/e30x3) (uf (OH /e:0x3) + 
[(1 — Pr)/Pr](oh €30%3)} — Rv;H') (le) 
If the outer flow is irrotational (or practically so), 
and if the lines (x. = constant, x3; = 0) are chosen coin- 
cident with the streamlines at the outer edge of the 
layer, the boundary conditions are 


X3 => 0 V; = V> = V3 = V1V3 = V2V3 = pv3 = 


H = H, 
Vi> Vy; Ve m0; Ho H,; 


x37 @© 
> 


—>0; pv; —> 0; | 


v3h1’ —> 0 | 


0103 —> O; V2V3 


Also, the momentum equations governing the inviscid 
flow become 


Vi-(OVie €10X}) = —(J R,)(OP, €;0%}) (3a) 


(1 R,.)(OP, €20X2) (3b) 


(Vi? /e1€2) (Oe1/ 0X2) 


The determination of the mean motion from Eqs. 
(la) through (le) and the boundary conditions (2) re- 
quires the a priori knowledge of relations connecting 
the turbulent stresses and the distributions of mean 
flow properties. Relationships analogous to those 
valid in laminar flow [s = u(def V), K = k grad T | were 
developed in the past on the basis of phenomenological 
theories, which regarded the statistical effect of turbu- 
lence on the mean flow as similar to that of viscosity. 
The accuracy of these theories has recently been ques- 
tioned by many authors. In contrast to the phenom- 
enological point of view, it is now generally accepted 
that the relationships in question should be based on a 
priori determination of the factors governing the motion 
in the boundary layer; for example, Townsend* has 
obtained one such result by assuming that the incom- 
pressible turbulent boundary-layer motion on a flat 
plate is controlled by a group of large eddies similar to 
those encountered in wake flow. In any case, rela- 
tions between the turbulent transports and the dis- 
tributions of mean flow properties can be written im- 
plicitly, without a specific knowledge of the governing 
factors; this point of view is adopted in the following 
analysis of the system of equations of motion. 

In accord with the above considerations we let 
formally (7,7 = 1, 2, 3) 


vo, = v0,(V,, P, H) l 
pt’s = pv3( Vi, rot Hf) (4) 
H’v; = H'v,(V, P, m' 


Thus, a complete (though implicit) system of equa- 
tions is available for the problem on hand; the relevant 
equations are Eqs. (la)—(le) and (4). 

If we denote by 


x = (1 €1€2) (Oey Ox2) 


the component of the outer streamline curvature in the 
plane tangent to the surface of the body, and if we 
combine Eqs. (lc) and (3b), the momentum equation 
in the x, direction can be rewritten as 


Vi(OV2/e;0x1) + V2(OV2/e20x2) + 
[V3 + (pv3/R) ](OV2/e30x3) + (Vi V2/ee2) X 
(Oe2/O0x;) — (1/R)(O0/e30x3) [u(OV2/e30x3) — 
Rovw3) = Vi.2«[(Vi2/Vi.2) — (R,./R)] (5) 


For either « ~ 0 or [(Vi? /Vi.2) — (R./R)| > 0, the 
right-hand side of this equation is equal to zero. The 
former condition is satisfied by quasi-two-dimensional 
flows; the latter condition is closely approximated 
when the wall is highly cooled and V;,?/(/7, — Hy) <1, 
that is, for the category of hypersonic flows considered 
in this paper.f| Without inquiring about the details of 
the relation 

VoV3 = v:03(V;, P, AH) (6) 


¢ 


we then distinguish between alternatives; (a) the 
right-hand side of Eq. (6) is homogeneous in Vz; (0b) 
the right-hand side of Eq. (6) is not homogeneous in V2. 

In the case (a), and under the subject hypersonic 
conditions, the x. momentum equation is homogeneous 
in V2; in view of the homogeneous boundary conditions 
(2) we may then conclude that V2 = vw; = 0. Thus, 
the problem becomes two-dimensional in the stream- 
line coordinates, and, in principle, can be solved by 
using the remaining equations in the system. 

In the case (b), it is convenient to consider the system 
of Eqs. (1a), (1d), (le), (4), (5), and the equation 


V1(0/e,;0%1)(Vi V2) + V2(0/e20x2)(ViVe) + 
[Vs + (pv3/R) ](0/e30x3)(Vi Ve) + (Vi? Ve/ere:) &X 
(Oe;/Ox2) + (V2?/e:e2) (Oe2/0%1)(Vi — V2) = 
—(V2/R)(OP/e,0x1) + ViVie? «[(Vi?/ Vi-2) — 
(R,/R)|] + (Vi/R)(0/e:0x3) X 
[(u/e3) (OV. /Ox3) — Rows) + 
(V2/R)(0/e30x3) [(u/es)(OVi1/Ox3) — Ro,v3|) (7) 


The latter results when the x; momentum equation is 
multiplied by V2 and added to the x, momentum equa- 
tion multiplied by V;. Under the hypersonic condi- 
tions considered here, Eqs. (5) and (7) are homogeneous 
in V2 and 7v;3, and with the homogeneous boundary 
conditions (2) admit of the solution V2 = v2; = 0. The 
ensuing two-dimensional problem may, in principle, 
be solved by using the remaining equations in the sys- 
tem mentioned above, the x, momentum equation being 
replaced by the nonhomogeneous equation 


+ This statement is verified in the section on boundary-layer 


flows depending on two space variables. 
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Laminar heat transfer on a blunted cone at Ma = 6 
Distribution on the windward side 


Fic. 1. 
and angle of attack a = 15°. 
of the pitch plane (¢ = 0°). 


V23 (Vi, F, ff) = 0 


Thus, we have shown that in the presence of highly 
cooled surfaces and of moderate Mach Numbers of the 
outer stream, the crossflow and the pertaining Rey- 
nolds stresses are negligible even for large transverse 
pressure gradients. It may be of interest to notice 
that the same result could readily have been obtained 
by assuming that the components of the mean velocity 
and of the Reynolds stress parallel to the surface of the 
body have the same direction—..e., 


V2 V; = Vov3 ‘D103 (8) 


Under hypersonic conditions, the x, momentum equa- 
tion [Eq. (5)] is then homogeneous in V2, and the 
desired result follows. 

In conclusion, the turbulent boundary layer on the 
nose portion of an arbitrary three-dimensional blunt 
body in hypersonic flight can be studied by two-dimen- 
sional methods used in connection with a streamline 
coordinate system. In the absence of accurate meas- 
urements of pressure and turbulent heat-transfer dis- 
tributions on such bodies, we cannot presently assess 
the range of applicability of the proposed approxima- 
tion. However, for either zero or favorable streamwise 
pressure gradients,} an indirect estimate to this effect 
may be extrapolated from available comparisons be- 
tween theory and experiment in laminar flow; this will 
be done here. 

For the specified category of pressure gradients it is 
known that, due to the larger shearing stresses, smaller 
three-dimensional effects can be expected for turbulent 
layers as compared to laminar layers subject to the 
same boundary conditions; thus, the hypersonic ap- 
proximation should yield more accurate results in the 
turbulent case than in the laminar case. A simple 
and rapid method for analyzing laminar boundary 
layers in the presence of highly cooled surfaces and of 


7 These are the pressure gradients to be encountered on prac- 
cal configurations for hypersonic flight. 
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moderate Mach Numbers of the outer stream was pre- 
viously described in reference 9, where a comparison 
with some exact solutions was also presented. More 
recently, a comparison with experimental data ob- 
tained at the Polytechnic Institute of Brooklyn" has 
been carried out. The configuration considered was 
a 20° half angle blunted cone with a spherical nose, 
tested at = 6, a = 15°, H,/H, ~ 0.5. Satisfactory 
agreement has been obtained between theory and ex- 
periment as shown in Fig. 1 for the streamline on the 
windward side of the pitch plane, and in Fig. 2 for the 
streamline inclined 15° to the pitch plane at the stag- 
nation point. These results indicate that the hyper- 
sonic approximation may have a rather wide range of 
applicability. 

In the following section a method of analysis ap- 
propriate to the subject category of turbulent bound- 
ary-layer flows is described. 


(3) Method of Analysis 


Inspection of the Eqs. (la)—(le) indicates that, when 
the crossflow and the pertaining Reynolds stresses are 
negligible, the only remaining parameter representative 
of the three-dimensional effects is the length element 
in the direction normal to the streamline, which appears 
in the continuity equation. Upon determination of 
this geometric parameter from inviscid flow consid- 
erations as described in reference 9, the evaluation of 
the boundary-layer characteristics can proceed inde- 
pendently along any selected number of streamlines. 
In this stripwise analysis one can let 


and, thus, rewrite the equations governing the steady 
mean motion in the form 
(oO Ox) [eoRV1] -- (Oo Ox3) [eo( RV 4 pvs) | = 
Vi(OV OX) a [V3 of (pv; R)|(0V, O23) oS 
—(1/R)(OP/Ox) + (1/R)(0/Ox3) X 
[u(OV1/Ox3) — Rvs] 
P = P, + Of6] 
V\(OlT/Ox1) + [V3 + (pv3/R)](OH/Ox3) = 
(1/R)(0/dx3)(uf (OH /dx3) + 
[1 — Pr)/Pr](Oh/dx3)} — R 2311’) 


$ (9) 


Eqs. (9) are formally identical to those encountered in 
the analysis of turbulent boundary-layer flows depend- 
ing on two space variables; for example, substitution 
of the radial coordinate 7 in place of e: yields the form 
appropriate to axisymmetric problems. The same is 
true for the boundary conditions; if the vorticity in the 
outer flow is neglected, these are 


«3 = 0 V1 = V3 = 0; fn = H,,) 
H —> H, \ 


(10) 


x3 > © Vs > Vas; 


If we assume that the same flow mechanism holds 
locally on an axisymmetric body as on a two-dimen- 
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sional body, we may extend to the subject problem the 
correlation between compressible and incompressible 
turbulent boundary layers recently proposed by Mager.*® 
The key to Mager’s transformation is the assumption 
that the turbulent process is unaltered by the com- 
pressibility effects due to internal heating of the fluid. 
Thus, one postulates that (a) at any point in the flow 
field the turbulent stress is a property of the local mass 
element of fluid, and (b) the turbulent shear associated 
with a mass element is invariant under a transformation 
which leaves the stream function unchanged. 


The choice of the reference conditions is particularly 
important in establishing the correlation between com- 
pressible and incompressible turbulent flows. If the 
effect of the internal heating on the stress transforma- 
tion is to be accounted for correctly, the stagnation en- 
thalpy distribution across the boundary layer must be 
conserved. For an arbitrary pressure gradient, this 
requirement leads to choosing the reference conditions 
at the stagnation value, when Pr = 1 and the wall is 
adiabatic. On the other hand, when heat is transferred 
at the wall, the reduction of the equations to the in- 
compressible form is only possible for zero pressure gra- 
dient and uniform wall temperature. The latter con- 
clusion, first reached empirically by several authors 
(e.g., Eckert!®), is corroborated by the following con- 
siderations: If, according to the hypotheses, the 
effects of internal heating do not alter the quality of the 
turbulent process, we can still view the compressible 
turbulent layer as a combination of two mechanisms 
(the inner and the outer portion, respectively) con- 
nected by one primary element (the wall shear). Under 
this assumption the effect of the wall shear on the pro- 
file shape in the presence of uniformly smooth surfaces, 
zero pressure gradient, sufficiently high Reynolds 
Numbers (such that the wall stress gradient is negli- 
gible), and uniform wall temperature, must still be 
represented accurately by the correlation (Vi; — Vi,) + 
V;, vs. (3/5), Vig = (Tu / Pu)? being a friction velocity 
and £; being a suitably transformed coordinate along 
the normal to the wall.| Thus, the compressible layer 
must be self-preserving, and the pertaining distributions 
of density and enthalpy can be represented by constant 
average values. 


In the present paper we postulate that, in the pres- 
ence of highly cooled surfaces and of moderate Mach 
Numbers of the outer stream, the compressible bound- 
ary layer may be correlated to an equivalent incom- 
pressible one, even though strong pressure gradients 
and large heat-transfer rates may be encountered. 
After performing the transformation we shall show that 
the effect of the pressure gradient is negligible and, 
therefore, that the postulation in question is correct. 
With this objective in mind, and for /7,,/H, = constant 
over the entire body, we introduce in Eqs. (9) new inde- 
pendent variables defined by the expressions 

} The criteria for selecting equivalent quantities are discussed 


at greater length in Section (4). 
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¢ [ Vie Rete 1 ‘. Vi, { R ] (11 
~ = F /2 7 Xx ; $3 = 9 X3 y, 
”” Innit HH’? Jo Ry 


where the subscript 0 identifies thermodynamic proper- 
ties in the reference state, to be selected consistently 
with the ratio H7,,/H,. Also we set forth transformed 
mean velocity components V,* and V;* and trans- 
formed Reynolds stresses Ry 2\7;* in accordance with 
the aforementioned conditions that the stream function 
y, defined byt 
Roro(Ow Ox) = —e(RV; + pus) ] 
> (12) 
Roro(Ow Ox3) = €o RV, \ 


and the turbulent shear associated with a mass element 


be invariant under the transformation. Thus, 
Vi* = (% €2) (Ow 0&3) = (H,' * Vie) Vi | 
V3* = —(ro/e2)(OW/dE1) = (R/Ro) X | 
{ (H.'*/ Vie) (Rewo/Rete) (Vs + (pv3/R)] — + (18) 


V(Ox3 0£;)} 
Rowvv3* = (77, Vie") (Rome Ru.) R V3 


Finally, we consider the energy equation in terms of the 
enthalpy ratio 


H* = (#7 — H,)/(A. — Ae) (14a) 


and of the pertaining transformed turbulent transport 


Ry 23H’* = [1/(H. — Hy) \(H2/?7/Vi.) X 


(Rouo/ Rue) RozfT’ (14b) 


It should be noticed that the transformation in ques- 
tion is not unique when V;, = 0; hence, the method of 
analysis can only be applied away from a stagnation 
point, appropriate initial conditions being prescribed 
in the neighborhood of such a point. 

With the new dependent and independent variables, 
Eqs. (9) and boundary conditions (10) become, respec- 
tively, 


t roisareference length. 
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(0/0&1) (e2.Vi*) + (0/0&s)(e2V3*) = O 
Vi*(0V;* 0£;) + V3*(0Vi.* O£3) = 
(H./ Vie) (OVie/O81) [((R./R) — (V1?/ Vi”) ] + 
(0/dEs) [(Ru/ Rete) ¥o(OVi*/OEs) — v105*] 


> (15) 
Vi*(OH*/d8:) + Vs*(OH*/0Es) = | 
(0/2Es)((Ru/Rette) vol (OH*/Es) + | 
(1 — Pr)/Pr]{1/(He — Hy) \(Oh/dEs)} — | 
voHl’*) | 
&=0 Vit=Vi* = H*=0) 
; ; (16) 
&—> © V\* —> H,'?; H* — 1\ 


In the presence of highly cooled surfaces and of moder- 
ate Mach Numbers of the outer stream 


R./R ~ H/H, ~ H* (17) 


Under these conditions, and for Pr = 1, a solution of 
the energy equation is found in accordance with Rey- 
nolds analogy 

H/H, ~ (H — Hy)/(He — Hw) = Vi/Vse) 

: » (18) 

vil’ /(H, — Hy) = r1s/ Vi. | 
Indeed, for this solution the pressure gradient term in 
the momentum equation becomes vanishingly small. 
Justification is thus obtained for the proposition set 
forth in the three-dimensional analysis, namely, 
[(R,/R) — (Vi?/V1,”) | ~ 0, and for the postulation in- 
cluded in the present two-dimensional analysis, namely, 
that the effect of the pressure gradient is negligible for 
the category of problems considered. 

The reduction of the equations to incompressible 
form is completed by assuming a linear dependence of 
the viscosity on the temperature so that, at any station 
on the body, (Ru)/(R.u.) = 1 across the boundary 
layer. The following system of equations is then 


obtained: 
(0 £1) (€2 Vi*) + (O 0&3) (2 V3*) = 0 


Vi*(0Vi*/08:) + V3*(0OVi*/0E3) = (19) 
(Oo O€3) [vo(0V1* OF) — 203*] 


which is typical of turbulent boundary layers with zero 
pressure gradient in incompressible flows depending 
on two space variables. . 

It may be of interest to notice that a further reduc- 
tion of Eqs. (19) to a strictly two-dimensional form, 
by means of a Mangler type transformation, cannot be 
performed if the turbulent shear associated with each 
fluid element is to remain invariant. 

The solution of Eqs. (19) is obtained here by following 
closely the classical semiempirical methods generally 
used for strictly two-dimensional flows. This approach 
extrapolates to the problem in question the experi- 
mentally verified approximation that the same flow 
mechanism holds locally on axisymmetric bodies as on 
two-dimensional bodies. Thus, we write for the ve- 
locity distribution throughout the boundary layer 


Vi*/V;,* = (1/K) {In (F.," &3 vo) + 2(és 5) | (20a) 
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or V,* — Vi." —_ Vi, *f (és 6) (20b) 
as first proposed by von Karman.” Correspondingly, 
we obtain for the momentum thickness of the trans- 
formed boundary layer 


o* = o1(Vi,* 'Vi.*)6 “2 2(V;,* Vi") 76 (21a) 


with 


1 1 
y= [16 5)d(&3/6); C2 = f [f(&2/6) |?d (3/6) (21b) 
0 0 


Alternatively, in terms of the variable z and of the 
quantity D defined by 


z= K(Vi.*/Vi,*); In D = f(1) (21c) 
we obtain 
o* = (v/DV,,*) [cq — (Keo/z) Je’ (21d) 


Substitution of the expression (21d) in the integral 
momentum equation 


(1 /e2) (d/d&,) (e29*) = 7,,*/RoVi-*? (22) 


and subsequent integration? yield the following law for 
the transformed skin-friction coefficient c,*: 


In c,* + In | Vie*/2c,%) (1/e2) f eats = 
V/2Kc/* 7“ (9%) 


where £1. designates the station at which z = 0 (e,> 0). 
Incidentally, this law includes Van Driest’s result! 
that, for a turbulent boundary layer, the cone skin- 
friction coefficient is the flat plate coefficient evaluated 
at one half the Reynolds Number on the cone.  Fi- 
nally, from Eq. (13), which establishes the turbulent 
stress transformation, we obtain for the actual skin- 
friction coefficient 


Cr = (ue/po)es* (24a) 
and for the actual heat-transfer rate 
Qo = Pr-*!® (H, — Hy) (ue/uo)(ReVie/2)e;* (24b) 


Eq. (24b) is usually recast in terms of Nusselt and 
Reynolds Numbers (Nu and Re) referred to either local 
or stagnation conditions of the outer flow; here we 
choose the second alternative—namely, 

Nu = [Quw/(He — Hy) \ro(cp/k) es} 


» (2) 


Re = H,'*7)(R, Me)s(P RH.) s' | 


so that for each configuration only one calculation is 
necessary within the range of Mach Numbers for 
which the pressure distribution is insensitive to the free- 
stream Mach Number. Naturally, such as extrapola- 
tion does not include the effects of nonunity Lewis 
Number” and of deviations from thermodynamic 
equilibrium. 

Since it has been shown that the equations governing 
the problem in question are reducible to those pertain- 


7 Only the highest order terms in z are retained for z > 1. 
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ing to a flat plate, the appropriate nondimensional 
heat-transfer parameter is (Nu/Re”*). The results 
presented by Libby and Cresci'* '* and by Rose, Prob- 
stein, and Adams” indicate that turbulent heat-trans- 
fer measurements on blunt-nosed bodies over a wide 
range of enthalpy ratios are accurately correlated by 
Nu ~ Re*®. Therefore, the heat-transfer equation is 
recast into the following final form appropriate to a 
highly cooled turbulent boundary layer subject to a 
pressure gradient: 


(1/2)(Nu/Re**)!"? In | arr'"\Ke* X 


| 


Ci," “ Vie) (Res Re) (Mes Me) (To €2) X 


f (V1. /H-'!?) (Re/ Res) (Me/ Mo) (€2/%e) (dx1 r) | aa 


(Nu/Re*®)*'? In [((Nu/Re*®)!?] = 
[Pr!3 Re! ®(K?/4)(Vie/F eo!) (Re/ Res) (oMe/ oo) X 
(R.H,/P),''?]'"* (26) 
with K = 0.40, D/aq = 5.24 
as determined by various authors on the basis of Eq. 
(23) and of available experimental evidence for incom- 
pressible flows. 
In the following section the accuracy of the method 
proposed here is assessed by comparison with experi- 
mental data and with predictions of cther known 


methods. 


(4) Comparison With Experiments 


All the analyses of turbulent boundary layers involve 
postulations confirmed either a priori or a posteriori by 
agreement with experiment. For the hypersonic case 
considered in the present paper, several extensions of 
classical low-speed methods have been presented re- 
cently.'* In all these analyses assumptions are re- 
quired for the skin-friction law, for the form factor, 
and for the skin-friction heat-transfer relation. Which 
of the assumptions is really crucial cannot be determined 
a priori; neither can a comparison with experiments, on 
the overall basis of predicted versus measured heat- 
transfer rates, be very illuminating in this respect, due 
to the limited amount of data available. 

The method presented here is based on the one crucial 
choice of the appropriate reference conditions to be used 
in the correlation of the turbulent stresses between in- 
compressible and compressible flow on a flat plate in 
the presence of heat transfer; thus, a distinct advantage 
is obtained in that the sole new empirical parameter, 
namely the aforementioned reference state of the fluid, 
can be determined, once and for all, by comparing pre- 
dicted versus measured heat-transfer rates on flat 
plates, or on equivalent configurations, for various 
values of H7,,/H,. 

An exact rule for the determination of the reference 
state probably cannot be derived on the basis of the 
experimental information presently available. How- 
ever, inspection of the results (in particular, the corre- 


lation Nu ~ Re**) presented in references 13, 14 and 
15, for wall enthalpy ratios /7,/H, < 0.70 and outer 
flow Mach Numbers /, < 2.5, suggests the following 
comments: 

(a) Within the range of the measurements the refer- 
ence state depends very slightly, if at all, on the wall 
enthalpy ratio. 

(b) A choice of the reference conditions based on 
Eckert’s reference enthalpy method leads to predic- 
tions higher than the measured values. 

These remarks are in accord with the measurements 
of Lobb et al.,'* who observed a pronounced upward 
shift of the turbulent portion of the boundary-layer 
profiles with increasing heat transfer. Under these 
conditions, the turbulent shear is important in a re- 
gion where the stagnation enthalpy is approximately 
constant; thus, a transformation based on the stagna- 
tion conditions of the outer stream should still be ap- 
propriate in accord with comment (a). 

Some discretion must be used in carrying our con- 
clusions further. If the stagnation conditions of the 
outer flow represent the appropriate reference state for 
wall enthalpy ratios /H/,,/H7, < 0.70, the relation be- 
tween actual and transformed local skin-friction co- 
efficients [Eq. (24a)] becomes independent of the wall 
enthalpy ratio. This result is in contrast with physical 
intuition and with experimental evidence, which indi- 
cate significant increases in skin friction when the wall 
temperature decreases, other flow conditions remaining 
unchanged. The anomaly may be interpreted by 
correlating the wall stresses in terms of the slope of the 
mean velocity profile at the wall to obtain the equation 


i (Ret Rete) (Me Mo) c,* 


which reduces to (24a) only when R,u, = Ru, = con- 
stant across the boundary layer, as assumed in the 
transformation. It is well known that, particularly 
for large degrees of cooling, R,,u, is significantly different 
from R,u,; therefore, skin-friction predictions based on 
Eq. (24a) are inaccurate. Exactly the same situation 
arises in the analysis of laminar boundary layers. 
Also in this case the assumption Ru = R,u, = constant 
leads to the aforementioned difficulty for skin-friction 
predictions; however, satisfactory heat-transfer esti- 
mates are obtained thereby. Thus, while the assump- 
tion Ru = constant may be retained in the analysis of 
laminar flows, different values of the constant must be 
used in the final equations for skin friction and heat 
Criteria for selecting these reference values 


21 


transfer. 
have been presented by Fay and Riddell. 

The heat-transfer measurements (to be discussed 
further below) reported in references 13, 14, and 15, 
and the skin-friction data presented in reference 207 all 
indicate that, for identical boundary conditions, the 
same values of Ru = constant apply for turbulent as 


+ Sommer and Short applied the 7’ method of Rubesin and 
Johnson to correlate turbulent skin-friction data obtained in the 
presence of negligible pressure gradient and severe heat transfer 


at the wall. 
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Fic. 3. Experimental variation of Nusselt Number with 
Reynolds Number at a given station (s = 0.69) on a spherical 
body (from reference 13). 


for laminar flow, in the presence of even moderate wall 
cooling. This evidence corroborates the conjecture 
set forth above that only minor changes of flow proper- 
ties are encountered in the outer portion of a turbulent 
boundary layer under the subject conditions; indeed, 
significant changes of density and enthalpy must be 
limited to the sublayer if laminar criteria are to apply. 
Since the choice of the reference state is essentially un- 
important for laminar flow, we conclude that a corre- 
lation based on the outer stream stagnation condition 
is generally appropriate for the range of enthalpy ratios 
considered. However, the transformation (11) should 
be generalized to the form 


f= | (Vie/He!!?) (Retr/Rowo) dy 


x10 


& = (Vie TT ,}"*) { (R Ro) dx3 
( 


) 


(27) 


with R,u, depending on the wall enthalpy ratio, if skin- 
friction estimates are of interest. The corresponding 
transformation of the local skin-friction coefficients is 


Cy = (Ryu, ’Rewo)c s* (28) 


In the light of this discussion, it appears that all the 
analyses of hypersonic boundary layers based on ex- 
tensions of classical low-speed methods cannot supply 
accurate results unless a suitable stretching of coordi- 
nates, which accounts for the compressibility effects, 
is included in the computation of the local Reynolds 
Number. The stretching is provided by the change 
of variables suggested in Section (3). This point of 
view is corroborated by the following comparison with 
experiments. 

In the absence of accurate pressure and turbulent 
heat-transfer data for blunt bodies in three-dimen- 
sional hypersonic flow, the comparison is limited to 
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three axisymmetric configurations on which measure- 
ments were made by Libby and Cresci.'* 14 The 
models were formed from part of a sphere and from a 
sphere-cone combination tested in air by the shroud 
technique.'* No special pressure distribution was 
imposed and, therefore, the results can be considered 
to be of a general nature. Ratios of stagnation en- 
thalpy to surface enthalpy varied from 2.0 to 1.4 dur- 
ing the testing period, since a transient technique of 
heat-transfer determination was employed. A typical 
set of measurements at a given station is shown in 
Fig. 3. In Figs. 4, 5, and 6 the fully developed heat- 
transfer data are presented together with the pertain- 
ing pressure distributions and with the predictions ob- 
tained by several methods as listed below. 

(1) The method proposed by Rose, Probstein, and 
Adams in reference 15. 

(II) A method based on flat-plate skin-friction law, 
Reynolds analogy, and form factor 17, = —1 (ef. 
Libby and Cresci'’). 

(III) The method proposed here. 

Methods (I) and (II) are based on similar extensions 
of classical low-speed analyses, the only difference being 
in the constant values assumed for the form factor 
(H,. = Oin reference 15). The pertaining heat-transfer 
equations are no simpler than the equation suggested 
here [Eq. (26)]. The calculations have been carried 
out by assuming a turbulent stress transformation 
based on the outer stream stagnation conditions, and 
by setting z = 0 in the neighborhood of the stagnation 
point. The effects of vibrational excitation of the air 
at high temperature have been considered in evalu- 
ating the local flow properties. 

Fair agreement with experiments and improvements 
on the other predictions are obtained by the method 
of this paper. Even though the comparison is acknowl- 
edgedly limited as to the quantity and accuracy (prob- 
ably +10 per cent) of the data, the subject approach 
to the problem of turbulent heat transfer appears to be 
reasonably substantiated. The present comparison 
with data obtained under conditions of moderate wall 
cooling represents an unfavorable test for the pro- 
posed method of analysis, which was developed for 


conditions of high cooling. The method appears to be 
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applicable for all wall enthalpy ratios equal to or lower 
than those considered here. These limits of validity 
approximately coincide with the limits of the “cold 
wall” approximation in laminar boundary-layer theory. 

Some comparative considerations are suggested by 
the discrepancy between the present predictions and 
those of method (I), which was also proposed for highly 
The two methods essentially exhibit 
(a) no stretching of coor- 


cooled layers. 
the following differences: 
dinates was introduced in reference 15 to account for 
the compressibility effects, and (b) the form factor was 
assumed /7, = 0 in reference 15, while it is 7, = —1 
in the present analysis. 

Inspection of Figs. 4 and 5 indicates that item (a) 
is the important one, the effect of the form factor being 
evidenced by the difference between the predictions of 
methods (I) and (II). The influence of the coordinate 
transformation is evidenced when the characteristics 
of the transformed incompressible layer are described 
by a 1/7 velocity profile and by Blasius skin-friction 
law, as proposed in reference 15; the following expres- 
sion is then obtained for the heat-transfer parameter 
in terms of Nu, and Ke, based on local external condi- 


tions: 


Nur “ l Mo de 
- = 0.029 Pri" - x 
Re,*” Viel emcee? * \ped X 


x1 ]-1/5 
[ Vi-R Mees” x, | (29) 
J0 


This equation differs from that proposed by Rose, 
Probstein, and Adams." The satisfactory comparison 
between theory and experiments reported by these 
authors can be attributed to the fact that their meas- 
urements on sphere-cylinder models were confined to 
an area one or more diameters back from the shoulder, 
in order to attain a zero pressure gradient flow; under 
these conditions Eq. (29) and the corresponding one 
in reference 15 yield approximately the same predic- 
tions, the differences in the analytical expressions being 
obliterated by the length of run at constant external 


conditions. 
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(5) Concluding Remarks 


Two-dimensional and three-dimensional boundary 
layers in the presence of highly cooled surfaces and of 
moderate Mach Numbers of the outer stream have 
been considered. These conditions are representative 
of those encountered on the nose portion of blunt bodies 
and in the neighborhood of the blunt leading edge on 
moderately sweptback wings. 

For problems depending on two space variables, a 
method of analysis has been established based on an 
extension of the compressible-incompressible turbulent 
boundary-layer transformation due to Mager. The 
selection of the appropriate reference state for this 
transformation has been discussed. On the basis of 
available experimental evidence it has been suggested 
that, in the presence of wall cooling, the outer turbu- 
lent portion of the boundary layer exhibits a very 
small variation of stagnation enthalpy and, therefore, 
that the stagnation conditions of the outer stream rep- 
resent appropriate reference values over wide ranges 
of wall enthalpy ratios. The accuracy of the proposed 
method of analysis has been assessed by comparison 
with experimental data; satisfactory agreement with 
the measurements and an improvement with respect 
to the predictions of other known methods have been 
obtained. 

For the three-dimensional problems it has been shown 
that, under the subject hypersonic conditions, the cross- 
flow and the pertaining Reynolds stresses are negligible 
even for large transverse pressure gradients. There- 
fore, the boundary-layer characteristics and, in par- 
ticular, the heat transfer at the wall may easily be 
estimated by using an orthogonal curvilinear coordi- 
nate system with the streamlines of the inviscid flow as 
one family of coordinate lines. Expressions for the 
evaluation of the heat transfer [Eq. (26)] have been 
provided depending on the geometrical characteristics 
of the coordinate system; these can be determined 
either from pressure measurements or from theoretical 
investigations of the inviscid flow field. The range of 
validity of the proposed cold wall approximation has 
been estimated qualitatively on the basis of recent 
measurements in laminar flow, the argument being 
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that, for either zero or favorable streamwise pressure 
gradients, the larger shearing stresses encountered in 
turbulent boundary layers, as compared to laminar 
layers, allow smaller three-dimensional effects. 
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Wings With Minimum Drag Due to Lift 
in Supersonic Flow 


(Continued from page 20) 


lift seems possible when compared with an untwisted, 
uncambered wing of the same planform. The greatest 
reduction in drag occurs for wings with both edges 
swept back and with a subsonic leading edge. A 
characteristic feature of these wings is the upturned nose 
of the root section. 
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Flexural Vibrations of a Thick-Walled 
Circular Cylinder According to the Exact 
Theory of Elasticity! 


JOSHUA E. GREENSPON* 


Summary 


This paper gives results of the exact theory for flexural or non- 
axially symmetric vibrations of a thick-walled cylindrical shell 
having freely supported ends. The results also apply to flexural 
wave propagation in infinitely long cylindrical shells. Curves of 
phase velocity and frequency as a function of wavelength ratio 
and wall thickness ratio are given. Included are some compari- 
sons between the exact theory, the Timoshenko beam theory, and 
shell theory. 


Symbols 

a = outside radius of cylinder 

ai = inside radius of cylinder 

l = length of cylinder 

p = mass density of cylinder material 

E = modulus of elasticity 

G = shear modulus 

v = Poisson’s ratio 

Ur = radial displacement 

ug = tangential displacement 

Uz = longitudinal displacement 

6 = cubical dilation 

Wr, WO, Wz = components of rotation 

C, = velocity of a dilatational wave 

C2 = velocity of a rotational wave 

r, 0, 2 = cylindrical coordinates 

rr, 70, 12, 

22, 20, 00 = eX ymponents of stress (Love’s notation) 
€rr, €r8, Erz, 
€22, €29, €99 = Components of strain (Love’s notation) 

A, ue = Lame’s elastic constants 

m = number of axial half-waves in vibration pattern 

n = number of circumferential waves in vibration 
pattern 

R\(r) = a function of r associated with the radial dis- 
placement 

Rr) = a function of r associated with the tangential 
displacement 

R;(r) = a function of 7 associated with the longitudinal 
displacement 

p = natural frequency 

a = a;/ao, thickness ratio 

B = mrdo/l = mdo/X, length ratio 

do = outside diameter 

rN = wavelength 
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ki’ = shear constant (such that Q = k;’AGy) 

Q = shear force 

A = cross-sectional area 

¥ = angle of shear at neutral axis 

¥ = p/(mn, DVG, p = C/C2ratio of phase velocity 


of wave to velocity of a pure rotational or 
shear wave; or ratio of frequency to first 
mode frequency of torsional vibrations 


VA = pary/ (1 — »2)/E = ¥8x/(1 — »)/2 


Introduction 


| gues 1 gives basic solutions for the flexural 
vibrations} of a thick-walled cylindrical shell hav- 
ing freely supported ends, along with results for long 
wavelengths. These solutions, which also represent 
traveling harmonic waves on an infinitely long hollow 
cylinder, are an extension of the solutions of references 
2, 3, and 4 for a solid cylinder. 

In this paper, the basic solutions of reference 1 are 
used to compute extensive phase velocity and frequency 
curves. Since the analysis is based on the three- 
dimensional theory of elasticity, it should serve a 
number of purposes, among which is the evaluation of 
the approximate shell theories. 


General Theory 


The equations of motion governing the free vibra- 
tions of an elastic body are the well-known three- 
dimensional equations of elasticity (see reference 5): 


t Also called nonaxially symmetric vibrations 





Fic. 1. Sketch of cylinder. 
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07u,/Ot? = C,2(06/Or) + 2C2? X 
[(Ow,/O0z) — (1/r)(Qw,/08) | 
d2n,/dt2 = (C,2/r)(08/00) + 2C2? X a) 


[(Ow,/Or) — (Ow,/O0z) | 
0°u,/Ot? = C\?(06/0z) + 2(C2?/r) X 
[(Ow,, 0#) — (0/0r) (rw9) | 


in which 
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E (1 — v) a 1 


C? = Fs > C2? = 
p (1+ v) (1 — 2p) p 2(1 + »v) 


where C; is the velocity of a dilatational wave, C2 is the 
velocity of a rotational wave, and 
6 = (1/r)[O(ru;)/or}] + (1/r)(Oug/00) + (Ou,/dz) 
the cubical dilation 
)[(1/r)(Ou,/00) — (Oug/dz)] 
( 


9 
2)[(Ou;/Oz) — (Ouz/Or)] 
2){ [d(rue)/dr] — (du,/oo)} | tation about 


or = the components 





(1 
we = (1 of average ro- (2) 
i; 


wo: = 
the three co- 


ordinate axes 


r is in the direction of increasing radius (r = 0 is the equation of the cylinder axis), z is in the direction of the axis, 
6 is in the direction of the tangent to the circles of constant radius, and u,, u», u, are the displacements of the element 


in the three coordinate directions. 


The components of stress and strain in the cylinder may be expressed in terms of the displacements by the following 


expressions: 


Stress Components 


- 


rr = 6 + Qu (Ou,/Or), 


| 


76 = w[(1/r)(Qu,/00) + (Ou,/dr) — (ue/r)], 


u[(Ou,/Oz) + (Ou,/Or) |, 


II 


'S 


zz = Xd + 2Qu(du,/O-) 
= u[(1/r)(Ou,/00) + (Ou»—/dz) | (3) 
06 = 6 + Qul[(1/r)(Om_/d0) + (u,/r)]' 


2g 
> 
| 


\ and uw are Lame’s constants, given in terms of / and v by the following expressions: 


\ = Ep/(1 — 2r)(1 + »); 


Strain Components 


€rr = Ou,/OFr, 


€, = (1/r)(Ou,/00) + (Ou,/Or) — (u,/r), 


€r2 = (Ou,/0z) + (Ou,/Or), 


With respect to the hollow circular cylinder, several 
special types of motion may be considered. In this 
paper we are concerned with the flexural vibrations of 
the cylinder which is subject to the following boundary 


conditions: 


22 =U =u=0 at 2 ps4 


t (6) 
rr = v6 = rz =~ VO at r=4,7 = A 


Reference 1 gives a complete discussion of the func- 
tions that describe the flexural modes of vibration of a 
thick-walled circular cylinder. Briefly these functions 


are 


u, = R, (r) sin (mxz/l) cos n 0 cos p ft] 
Ug = R> (r) sin (maz/l) sin n 6 cos pt (7) 
u, = R; (r) cos (mmz/l) cos n 6 cos pt 


The functions Rk), Rs, and R; have been determined 
such that Eqs. 7 form a solution to the three-dimen- 
sional equations of elasticity. The notation for the 
displacements is given in Fig. 1. 

It can be verified that Eqs. 7 also represent super- 
positions of two harmonic waves traveling in opposite 
directions along the axis of an infinitely long cylinder 
with phase velocities equal to p/(m7z/l). Thus, all the 
curves of this paper can be used to compute phase veloc- 
ities in an infinitely long hollow cylinder. These phase 


pw = E/2(1 + »v) (4) 
€z2 = OUu,/O2 | 
€x9 = (1/r)(Ou,/00) + (Ou9/0z) ? (5) 
€99 = (1/r)(Ou,/080) + (u,/r) { 


velocities then can be used to compute group velocities 
of propagating disturbances in thick-walled cylinders. 

The frequency determinant (see reference 1) is 
formed by first equating the stresses on the outside 
and inside cylindrical surfaces to zero. This gives a 
set of homogeneous equations from which the fre- 
quency determinant is developed. 


Results 

The roots of the frequency determinant for n = 1, 2, 
3 were computed for a series of thickness and length 
ratios. These results are given in Figs. 3 to 9. Fig. 3 
shows the first branch phase velocity curves for a hol- 
low cylinder undergoing beam-type vibrations corre- 
sponding tom = 1. All curves are asymptotic to y = 
0.93, the Rayleigh wave velocity. For 0.7 < a < 
0.98 a characteristic dip is seen in the curves; this dip 
gets steeper as the shell gets thinner. For longer waves 
in the region where 6 < 1, the y parameter is very in- 
sensitive to thickness. However, for the intermediate 
values of 8, the effect of thickness is more pronounced. 
Only when 6 becomes very large will the curves come 
together again. Fig. 4 shows the curves of Fig. 3 
plotted for 0 < 8B < 5. 
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Fic. 4. Phase velocity parameter as a function of length ratio for 
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Fic. 5. Comparison of exact theory with thin shell theory of ref- 
erence 6 and Timoshenko beam theory (” = 1, v = 0.3). 











FLEXURAL VIBRATIONS OF A THICK-WALLED 


CYLINDER 








Phase velocity parameter as a function of length ratio for 
= 0.3 (first root). 


) 


n = 2, 


v 











Phase velocity parameter as a function of length ratio for 
= 0.3 (first root). 


n=3 


»v 





-—— 


| 


deme 








\\ 














ett 
bee 














SENN 


Modified frequency parameter as a function of length 


ratio for n 


9 
2,¥ 


0.3 (first root). 








40 JOURNAL OF THE AERO 














\M\ 
\ \\ 


























Na 
Fic. 9. Modified frequency parameter as a function of length 


ratio for n = 3, v = 0.3 (first root). 


Fig. 5 gives a detailed comparison between the exact 
theory, Timoshenko beam theory, and the thin shell 
theory of reference 6. For a shell with a = 0.98, it is 
seen that the thin shell theory is very accurate, exhibit- 
ing very small errors even for rather short wavelengths. 
For a = 0.9, thin shell theory is good for values of 8 
up to about 8, corresponding to (Outside Diameter + 
Wavelength) equal to about 2'/.. For a = 0.9 it is 
seen that Timoshenko beam theory is inadequate and 
does not give the correct form of the curve for 6 > 1. 
Here, a shear constant of 0.5 was used; this value of 
shear constant gave excellent results for 0 < 6 < 1. 
For a = 0.8, thin shell theory becomes inadequate 
when 6> 5. For a = 0.5, thin shell theory gives the 
proper form of the curve but is very inaccurate even for 
small values of 8. On the other hand, Timoshenko 
beam theory, using a shear constant of 0.75, gives rather 
good results. For a = 0.01, the thin shell theory is 
completely inadequate and the Timoshenko beam 
theory with a shear constant of 0.9 (as recommended 
in the literature) is excellent. 

Fig. 6 gives the set of first branch phase velocity 
curves for the first lobar mode corresponding to ” = 2. 
For all thicknesses, the phase velocities approach the 
velocity of Rayleigh waves as 6 tends to infinity. 
Here there is also a characteristic dip in the phase 
velocity curves of the thinner shells. However, the 
very thick walled cylinders (corresponding to a = 
0.01) do not exhibit this dip. Contrary to the 7 = 1 
case, the phase velocity for m = 2 tends to infinity as 8 
tends to zero. Fig. 7 shows the phase velocity curves 
for n = 3. 

Since the phase velocity tends to infinity as 6 tends 
to zero, for nm = 2 and 3, another parameter, V/A, can 
be introduced to define the frequency for smaller values 
of 8. This parameter has been used by Arnold and 
Warburton.’ Figs. 8 and 9 give ~A asa function of 
8. Itis seen that the frequency curve has a flat portion 
for small 8. As the shell gets thicker this flat portion 
extends over a larger range of 8. This flat portion indi- 
cates that the frequency is independent of the wave- 
length in that region. 


SPACE 
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Further Shell Theory Comparisons 


In a recent letter, Dr. G. B. Warburton (coauthor of 
reference 6) told of some approximate calculations he 
has made form = 2and1<6< 3.3. These show that 
the thin shell theory of reference 6 gives first branch 
frequencies accurate to 10 per cent for 1 < 8 < 3.3 with 
a = 0.9, 0.8 and for 1 < 6B < 2 with a = 0.7. 

There are a number of existing approximate thick 
shell theories, some of which are summarized excel- 
lently by Cooper and Naghdi.’? Judging by some 
comparisons between the results given in reference 7 
and the exact theory of this paper, it seems that the 
theory described by System (I) equations in reference 
7 offers a great deal of promise. Some rough calcula- 
tions based on the first mode curves of reference 7 for 
a = 0.9, n = 1, 2, 3 indicate that the results predicted 
by System (I) equations compare very favorably with 
the three-dimensional theory. Even y values for large 
8 (up to 8 ~ 20) are predicted rather well by reference 
7. Other approximate thick shell theories such as the 
one contained in reference 8 have not been examined 
in detail as yet. 

Some other results for nm = 0, 4, 5 and higher branches 
of the dispersion curves as well as mode shapes and 
stress distributions have been computed recently and 
are contained in reference 9. Also Gazis'! has pub- 
lished an extensive set at higher branch curves for n = 1 
and 2. His theory is equivalent to the theory pre- 
sented in reference 1. The results of Gazis, together 
with those of the author, should form a rather complete 
picture of the behavior of thick-walled circular cylinders. 
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Effectiveness of Radiation as a Structural 
Cooling Technique for Hypersonic Vehicles 


ROGER A. ANDERSON* anp WILLIAM A. BROOKS, JR.** 
Langley Research Center, NASA 


Summary 


The effectiveness of radiation as a structural cooling technique 
is examined for vehicles which are subject to large variations in 
heating over their surface. An analysis is made of the more 
significant effects of internal heat transfer within idealized struc- 
tures, both with and without external insulation on areas of 
greatest aerodynamic heating intensity. Selective use of insu- 
lation leads to structural temperature reductions substantially 
greater than those obtained without insulation. The results are 
presented in nondimensional form, and numerical examples 


illustrate applications to hypersonic flight vehicles. 


Symbols 
d = diameter or thickness of leading edge, ft. 
E = modulus of elasticity, lb. /in.? 
i] = function of d/l 
F,-» = interchange factor for radiation emitted by area a and 
intercepted by area b 
k = thermal conductivity, B.t.u./ft.-hr.-°F. 
l = leading-edge chordwise dimension, ft. 


heating rate, B.t.u./ft.?-hr. 
t = thickness of insulation layer, ft. 


S 
I 


T = temperature, °R. (unless noted otherwise) 
w = weight per unit area, Ib./ft.? 
€ = material emissivity 
e = function of emissivities (defined in appendix ) 
p = density of material, lb. /ft.* 
o = Stefan-Boltzmann constant, 0.1714 X 10°% B.t.u. 
ft.2-hr.-(°R.)4 
Cy = yield stress, lb./in.? 
Subscripts 
e = equilibrium 
iL = lower surface of wing 
O = outer surface of insulation 
Ss = stagnation area of leading edge 
U = upper surface of wing 
W = closing web of leading edge 
Introduction 


N HYPERSONIC FLIGHT, aerodynamic heating rates 
I can be quite high, with large variations over the 
surface of a vehicle. In a vehicle developing lift, for 
example, the heating rates in stagnation areas are more 
severe than those on the windward side of the lifting 
surfaces which in turn are higher than those on the 
leeward side. The resulting temperatures introduce 
Presented at the Thermal Stress Session, IAS 27th Annual 
Meeting, N.Y., January 26-29, 1959. 

* Head, Airframe Components Branch, Structures Research 
Division. 

** Head, Structural Configurations Section, Structures Re- 


search Division. 


structural problems in controlling thermal stress and 
distortion as well as in choice of materials for con- 
struction. 

The very existence of large variations in surface heat- 
ing, however, presents a number of opportunities for 
reducing structural temperatures in highly heated 
areas while raising those in cooler areas, thereby allevi- 
ating the more critical design problems. The purpose 
of this paper is to examine the changes in structural 
temperatures that are possible with internal transfer 
of heat between adjacent areas and with selective use 
of external insulation when radiation is the only means 
utilized for dissipating heat to the surrounding atmos- 
phere. 

This purpose is carried out by a steady-state analysis 
of the heating input and internal heat transfer in critical 
structural areas of winged vehicles designed for flight at 
speeds up to and including orbital velocity. The basic 
results are presented in nondimensional form and illus- 
trated with numerical examples. Implications of the 
results in regard to structural design and material 
selection are discussed. 


Configurations Analyzed 


In the region of the wing leading edge of a hypersonic 
glide vehicle, large differences in heating intensity can 
occur on adjacent surface areas. This situation is de- 
picted in Fig. 1 where the relative intensity of heating 
on three surfaces of interest, the wing lower and upper 
surfaces and a stagnation surface, is indicated by sizes 
of the arrows. 

A first approximation to the structural temperature 
at a given location is provided by the temperature 7, 
at which the local convective heating input g to an in- 
ternally insulated nonconducting skin of zero heat ca- 
pacity is balanced by the local radiation output. Radi- 
ation equilibrium temperatures can be calculated with- 
out reference to structural details and are the maxi- 
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Idealized aerodynamic heating input to structural area 
of hypersonic glide vehicle. 
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Fic. 2. Configurations used for heat-transfer analysis. 


mum possible temperatures that structures can attain. 
In this paper, these temperatures will be used as refer- 
ence values with which structural temperatures will 
be compared. It may be verified'~* that values of 7; 
for the more highly heated areas lie in the range 2,000° 
to 4,000°F., which suggests a serious problem in mate- 
rial selection for load carrying structures. 

Temperatures attained by the structure are influ- 
enced by internal radiation and conduction of heat and 
can be greatly reduced by the presence of external in- 
sulation and/or a source of internal heat capacity. 
An investigation of the reduction in structural tem- 
peratures by these combined means generally requires 
reference to a specific vehicle structure and a specific 
heating time history. The transient heating calcula- 
tions involved are often laborious. However, com- 
parisons of transient and steady-state calculations 
have indicated that during the gliding flight of winged 
re-entry vehicles, the velocity and altitude change 
slowly enough to permit the use of steady-state 
analysis with satisfactory results. In other situations, 
steady-state analyses of the heat transfer in idealized 
structures provide good insight to the effect of internal 
radiation. : 

Several idealizations are employed for this purpose, 
as shown in Fig. 2. 

Case I in Fig. 2 represents a wing cross section in 
which significant internal heat transfer occurs between 
essentially parallel surfaces of equal area. If the two 
surfaces see each other, as is likely to be the case in 
outboard areas of the wing, the transfer will occur 
principally by internal radiation. Case II represents 
the localized region of the wing leading edge where 
heat is transferred from the stagnation area by radi- 
ation and conduction to the wing upper and lower 
surfaces. Inasmuch as the heat-transfer rate in this 
region depends on the relative area of the skin elements 
involved, a structural parameter //d is involved. 

When a reduction in the temperature of the most 
highly heated skin element occurs because of internal 
transfer of heat to cooler areas, a further reduction in 
temperature can always be obtained when this element 
is externally insulated to restrict the inward flow of 
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heat. This approach is illustrated in Cases III and 
IV of Fig. 2 where insulation is applied only to the most 
highly heated areas, leaving the areas of lesser heating 
uncovered to act as radiators for internally transferred 
heat. Insulation conductivity and thickness are intro- 
duced as additional variables in the analysis. 

The effectiveness of this approach to radiation cool- 
ing of a structure is dependent on the existence of sur- 
face areas which are subject to much lower rates of con- 
vective heating than the average for the vehicle and 
which are in the immediate vicinity of the most highly 
heated areas. This situation occurs on the wing of a 
glide vehicle which is programmed to fly at moderate 
positive angles of attack during the period of most in- 
tense heating. Insulation of all external surfaces will 
not lead to any significant reduction in structural tem- 
peratures in a steady-state or long-time heating situ 
ation. Heat would simply flow into the structure at a 
reduced rate until internal temperatures reach values 
sufficiently high to force heat out through the insulation 
barrier at the coolest surface of the wing at a rate fast 
enough to balance the flow of heat in at the hotter 
surfaces. 


Uninsulated Structure 


Solutions for the skin temperatures of the uninsulated 
structures shown in Fig. 2 are derived in the appendix 
for several specific cases of internal heat transfer. 
The resulting equations and a discussion thereof are 
presented below, first for the case of the wing surfaces 
and then for the leading edge. 


Case I. Wing Surfaces 


The temperature of each skin element is conveniently 
expressed as a fraction of the radiation equilibrium 
temperature of the most highly heated element. Fol- 
lowing this approach, the temperatures of the upper 
and lower skins for the case of radiant transfer of heat 
between surfaces of equal emissivity can be written 


Ty P «. L= } {1 + (3 = €) (du qt) | (4 — e)} 1/4 (la) 
T1/T.,. 1 = \[(Qu/qu) + 3 — €)/(4 - e)}} * (1b) 


If perfect heat transfer (solid wing of infinite conduc- 
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Fic. 3. Effect of internal heat transfer on wing skin temperatures, 
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tivity) existed, the temperatures of both surfaces would 


be 

T/T., 1 = ((1/2)[(1 + (qu/qu) (2) 

and if no transfer of heat took place the temperatures 
would be 

Tu/Te, 1 = (Quv/qr)"" (3a) 

T,/T.,1 = 1 


(3b) 


Results calculated from these equations have been 
plotted in Fig. 3. The lower group of curves gives the 
temperature of the upper skin for the three types of in- 
ternal heat transfer, and the upper group of curves ap- 
plies to the temperature of the lower skin. For a glide 
vehicle flying at positive angles of attack, the appli- 
cable range of upper to lower surface convective heating 
rate ratio gy/qg, is between 0 and 0.3. In this range, 
significant reductions in lower surface temperature 
occur with internal heat transfer by radiation, and sub- 
stantial increases in upper surface temperature result. 
The temperature difference between the surfaces is seen 
to be slightly influenced by material emissivity, and 
decreases with an increase in emissivity. 

An indication of the effect which conduction of heat 
between surfaces may have on their temperatures 
may be obtained by consideration of the curves in 
Fig. 3 for radiant heat transfer and perfect heat trans- 
fer. Fig. 3 shows that radiation alone accounts for 
more than half of the possible changes in upper and 
lower surface temperatures that can occur with internal 
heat transfer. Investigation of a further reduction 
in the temperature difference 7, — 7y due to coupling 
conduction with radiation shows that, for a wing of 
reasonable depth and the usual amount of internal 
connecting structure, conduction has a_ negligible 
effect on the surface temperatures when 7,,, is at 
the 2,000°F. level and has only a minor influence when 
T.. ,isat the 1,000°F. level. 

Note that material emissivity determines the tem- 
perature level in Fig. 3 in a direct way through the re- 
lationship 

T., 1 = (Qz/ce)*!4 (4) 
A value of emissivity, as well as of gz, must be assumed 


to calculate numerical values of either 7, or 7, from 
values of T7/T,, , and gy/gz. To show the importance 
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Fic. 4. Plot of radiation law for selected values of emissivity. 
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temperature 


of emissivity in this calculation the radiation law 
[Eq. (4) ] is plotted in Fig. 4 for selected values of «. A 
curve for a value of « = 2 is included to indicate the 
lower limit of equilibrium temperature possible for an 
uninsulated skin An effective emissivity 
value of 2 is achieved when black-body radiation to 


element. 


absolute zero occurs from both surfaces of a skin ele- 
ment, a result given by the curve for perfect internal 
heat transfer in Fig. 3. 


Case II. Leading Edge 


Two modes of heat transfer are considered for this 
internal radiation. alone, and the limiting case of 
The solution for the 


case: 
infinite material conductivity. 
stagnation area temperature when heat is transferred 
to the wing upper and lower surface by radiation may 


be written 


(Ts F 6, s)* = 
f(d/l) + [1 — f(d/) \ (qu + gz) /2¢s)] (5) 


where f(d/l) is a function of leading-edge geometry 
given in the appendix. Similarly, for perfect con 


duction of heat away from the stagnation area, the 


solution is 
( Ts ) - l + 2(//d) [(qu + gz) /2qs]| 
a [1 + 2(1/d)] 


In both Egs. (5) and (6) the stagnation area temper- 
ature is given as a function of the ratio of the average 
of the upper and lower surface heating rates to the stag- 
nation area heating rate. A likely range for this ratio 
is G to 0.2; accordingly, this range is used to present an 
evaluation of Eqs. (5) and (6) in Fig. 5. The results 
are presented only for values of the leading-edge struc- 
tural aspect ratio //d up to five, inasmuch as simplifying 
assumptions in the analysis lead to temperature pre- 
dictions of decreasing validity as / d increases. Con- 
vective heating inputs and temperatures were assumed 
constant over each of the surface areas and black-body 
radiation was assumed to occur between areas. 

The major reduction in stagnation area temper- 
ature, due to internal radiation of heat to the cooler 
surface areas, is seen to occur at small values of //d. 
The percentage reduction in temperature is comparable 


(6) 
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Effect of lower surface insulation | on wing structural 
temperatures. . 
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to that shown in Fig. 3 for the lower surface of the 


wing. 


The dashed curves in Fig. 5 indicate the potential for 
temperature reduction due to material heat conduc- 
tion. An investigation of the importance of heat 
conduction in specific leading-edge designs* shows 
that, with a nearly solid leading edye, material temper- 
atures in the stagnation area can be reduced to about 
80 per cent of the radiation equilibrium value. Again, 
only a small //d is required because of steep temper- 
ature gradients associated with a high rate of stagna- 
tion heating and a finite value of material conductivity. 
From the present results and those in reference 4, it may 
be concluded that the solid curves of Fig. 5 give a real- 
istic estimate of the temperature reduction in a lead- 
ing edge in the form of a thin metal shell, because a 
very substantial material thickness is required to ob- 
tain an important benefit from heat transfer by con- 
duction. 


Insulated Structure 


Analytical expressions for evaluating the effect on 
structural temperatures of an insulation layer on the 
wing lower surface and stagnation area are derived in 
the appendix. These expressions are based on con- 
duction of heat through the insulation layer and internal 
radiation of heat to the cooler structural surfaces. The 
resulting equations and a discussion of the results are 
presented below. 


Case III. Insulated Wing Surface 


In the presence of an insulation layer, structural tem- 
peratures are conveniently expressed as a fraction of the 
radiation equilibrium temperature at the insulated 
surface. The subscript O (for outer) will be used to 
denote temperatures at the surface of the insulation 
layer. The three temperature ratios of interest, 7y + 
Te, 0, Tx/Te, 0, and To/T., 0, may be obtained from 
simultaneous solution of the following three equations 
which are valid for emissivity values near unity (see 
appendix) : 


AERO/SPACE 
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Tu/Te, 0 = (1/(2)"*)[(T1/Te, 0)* + (Qu/gz)\""* — (7a) 
To x qo = [1 — (T, Te. o)* + (Ty Ze or (7b) 
gr see, (To T,, o) (T/T. (7c) 


RTe,0 = (T/T, 0)* — (Tu/Te, 0)* 
The nondimensional insulation parameter tg,/kT., 0 
can be expressed in the following alternate forms 


tgr/kT., o = wqr/RpT., o = waoel., o°/kp (8) 


The alternate forms are particularly useful in weight 
analysis, inasmuch as w is the weight of insulation per 
unit area and kp is a measure of the efficiency of a par- 
ticular insulation material. 

Some solutions of Eqs. 
The temperature of the wing lower surface 7',, which is 
of primary interest, is given for a number of values of 
gu/9x. The temperatures 7; and 7o are given for 
values of gy/gz equal to 0 and 0.3 to give an indication 
of their trend in relation to 7,;. For a zero value of 
the insulation parameter, the values of Ty, and 7, are 
those given in Fig. 3. As the insulation parameter 
approaches infinity, the temperatures approach the 
following values 


T1/Te, 0 = Ty 
To aa. o= 1 (9b) 


(7) are presented in Fig. 6 


(qu/qx) "4 (9a) 


These values are closely approached for a value of 
wgr/kpT., o = 100. 

A comparison of Figs. 3 and 6 reveals that the re- 
duction in the structural temperatures 7,, and 7, ob- 
tainable with insulation on the hotter surface is very 
substantial, particularly when gy/gz is small. Insu- 
lation also reduces the temperature difference 7, — Ty 
for any value of gy/qz. Numerical values of 7, asso- 
ciated with changes in insulation for a constant radi- 
ation equilibrium temperature of 2,000°F. and « = 
0.9 are included in a table presented in Fig. 7. 

Values of t/k equal to 0.5 and 1.0 given in Fig. 7 
correspond to the practical range of values for insula- 
tion panels suitable for the lower surface of a wing. 
Associated values of 7, approach the limiting tem- 
peratures for a perfect insulation on the wing lower 
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surface. An evaluation of the insulation weight re- 
quired to obtain these temperature reductions is given 
in a subsequent section entitled “Weight Consider- 


ations.”’ 


Case IV. Insulated Leading Edge 


As noted previously, the internal transfer of heat by 
radiation from the stagnation area to the aft skin ele- 
ments of a leading edge of proper aspect ratio (//d = 2) 
takes place at about the same rate as transfer of heat 
from the lower to the upper surface of the wing. It 
follows, then, that the effect of an insulation layer in 
reducing structural temperature in the stagnation area 
can be estimated from Fig. 6. When Fig. 6 is used for 
this purpose, the temperature 7), is interpreted as 7's 
and the heating rate ratio gy gq, is taken to be (qy + 
gt) /2qs. The curves for Ty give an approximate indi- 
cation of the temperatures in the skin elements behind 
the stagnation area. 

Of interest is the fact that the refractory materials 
required for insulation of the stagnation area can be 
of much lower insulating efficiency (higher values of 
kp) and still produce percentage reductions in temper- 
ature equal to that achieved on the wing lower surface 
with an equal unit area weight of insulation. For 
example, if the radiation equilibrium temperature at 
the stagnation area is twice that at a point on the wing 
lower surface, a constant value of woe7, 9°/kp is 
maintained with an eightfold increase in kp. Because 
the convective heating rate ratio (¢y + @z)/2@s is also 
likely to be less than the ratio gy/gz, an even larger 
increase in kp can be tolerated and still maintain the 
same percentage drop in temperature through the 
stagnation area insulation. 

The table in Fig. 7 contains structural temperature 
values appropriate to the stagnation area of a re-entry 
glide vehicle. The temperatures given are for values 
of t/k that are one-tenth those used to calculate 
wing lower surface temperatures. In spite of this re- 
duction in insulating value, stagnation area tempera- 
tures approach their minimum possible values. 


Weight Considerations 


There can be at least two structural objectives in 
application of insulation to a vehicle. One might be 
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simply a reduction in the operating temperature level 
of the structure with attendant decreases in temper- 
ature differences between adjacent parts. A more 
comprehensive design objective would be an overall 
weight reduction. The successful achievement of either 
goal depends to a large extent on required insulation 
weight. 

The parameter kp is a useful index to the weight 
efficiency of a given insulation material. A com- 
parison of presently available materials and composite 
insulation panels on this basis is presented in Fig. 8 
for a range of temperature. Two distinct ranges of kp, 
differing by approximately two orders of magnitude, 
are indicated. The upper range is applicable to high 
density refractory ceramic oxides and contains those 
materials readily available for insulation of the stag- 
nation area. The lower range is occupied by low-den- 
sity fibrous insulation materials fabricated into various 
types of lightweight panels which can withstand the 
flow condition on a wing. The values of kp for these 
panels include the effects of structural attachments. 
Some of these panels have an exterior metal surface 
which limits the operating range of temperature. 

Using these lower values of kp, the unit area weight 
of insulation required to achieve desired reductions in 
wing temperature has been calculated and is presented 
in Fig. 9 A radiation equilibrium temperature of 
2,000°F. and a material emissivity of 0.9 have been 
assumed as nominal values for a hypersonic glide 
vehicle. The temperature of the lower skin of the wing 
associated with two extreme values of upper to lower 
surface heating rate ratio and a range of insulating 
material efficiency is given as a function of insulation 
weight. 

For a heating rate ratio of zero, the temperature of 
the lower skin may be reduced from a value of 1,770°F. 
without insulation to a value limited principally by a 
rapidly increasing weight of insulation. The impor- 
tance of a low value of &p for insulation of major areas, 
such as a wing, is readily apparent. For the higher 
value of the heating rate ratio, the maximum possible 
change in skin temperature is from 1,845°F. without 
insulation to 1,360°F. for perfect insulation. A struc- 
tural temperature of 1,400°F. can be achieved with an 
insulation weight of about 1 Ib. /ft.* of wing area. 











Qy oe 
a =O Teo 2000 °F 
2, €=09 
| kp 2 4 
\ 
\\ \ 
\ \ \ 
\ \ T#1,360 
INSULATION 
WEIGHT, | 
LB/SQ FT 
fe) 400 800 (200 (p00 2000 
LOWER SKIN TEMPERATURE, °F 
Fic. 9. Insulation weight required to reduce lower skin 


temperature. 





46 JOURNAL OF THE AERO 


Q| 


seed 


nm 
= 


1,500 2,00 


) 500 1,000 
TEMPERATURE, °F 


Effect of temperature on structural weight for various 
design criteria. 


Fic. 10. 


From Fig. 9, it may be concluded that wing struc- 
tural temperatures in the range from S00° to 1,400°F. 
are possible for a wing lower surface insulation weight 
of about 1 Ib./ft.* when radiation equilibrium temper- 
atures are at the 2,000°F. level. The actual temper- 
ature will depend upon control of the flight angle of 
attack (determining factor in gy /q,) of the vehicle. 

The question of offsetting the weight of insulation 
with a saving in the weight of the primary structure 
because of decreased operating temperature is one 
properly reserved for the designer of a particular ve- 
hicle. An indication of whether such a weight trade is 
possible, however, may be obtained on the basis of 
structural weight trends with increasing structural 
temperature. The simplest possible index to weight 
trends is the deterioration in structurally important 
properties of materials with temperature increase, pre- 
sented in Fig. 10. 

As can be seen from Fig. 10, structures designed solely 
on the basis of strength (yield stress, o,, important) 
undergo the greatest weight penalty with increased 
temperature, while structures designed on the basis of 
stiffness considerations (elastic modulus, /, important) 
undergo the least. The intermediate curve is appli- 
cable to structures in which failure due to compression 
buckling is an important criterion. 

The precise applicability of these weight trend curves 
to a hypersonic glide vehicle is beyond the scope of 
this paper. There is considerable evidence, however, 
that the design of such vehicles may be quite heavily 
influenced by aeroelastic considerations and by the 
availability and usability of minimum gages of sheet 
materials. In such circumstances, the lowest weight 
trend shown in Fig. 10 may be the most nearly appli- 
cable to the vehicle as a whole. 

Whatever the weight trend of an actual vehicle may 
be, it is clear that the severest test of the possibility 
for a favorable trade of insulation weight for struc- 
tural weight occurs with structures of very low unit 
area weight and designed by the use of criteria which 
lead to only small weight increases with increasing 


temperature. This type of structure has been assumed 


in preparing Fig. 11 where the combined weight of 
insulation and wing structure is shown as a function of 
the lower skin temperature. 


The external heating 
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conditions are the same as those given in Fig. 9; accord- 
ingly, the total weight at any value of lower skin tem 
perature is the sum of the insulation weights from Fig. 
9 and the wing structural weight. The wing structural 
weight is taken as 3 lbs./sq.ft. of plan-form area multi- 
plied by the variation of p/E with temperature given 
in Fig. 10. 

These calculations indicate that substantial struc 
tural temperature reductions may be obtained with no 
total weight increase over uninsulated structures, and 
that a significant reduction in temperature may be 
combined with a saving in total weight. A minimum 
weight and an optimum operating temperature occur 
for each value of insulation efficiency. These regions 
of minimum weight will shift toward lower temper 
atures in a vehicle of higher structural weight per unit 
area than that assumed in these calculations or in one 
whose structural weight increases more rapidly with 
increasing operating temperature. On the other hand, 
if minimum sheet thickness considerations dominate 
the structural design, only a very limited amount of 
structural material may be available for offsetting the 
weight of the insulation. In such a case, the use of 
insulation may be justifiable only on the grounds of 
reducing structural temperatures to the range of appli- 
cability of the commonly used nickel-base high-tem- 
perature alloys. 

The weight of insulating material required to protect 
the stagnation area of the structure against excessive 
temperature has been calculated for an assumed set of 
design conditions and is presented in Fig. 12. A range 
of kp, which varies from values typical of readily avail- 
able dense ceramics to values typical of more porous 
ceramics under development, is presented. 

Fig. 12 shows that the stagnation area unit weight 
of insulation is higher than the unit weight for the 
wing. However, the inside face or structural temper- 
ature reductions justify a high unit area weight over the 
relatively small total area involved in a leading edge. 
Temperatures can be reduced to values where a thin 
shell of a metallic material such as molybdenum (suit- 
ably coated for oxidation protection) is practical for 
the structure of the leading edge. The weight of such 
a shell with a ceramic nose is likely to compare favor- 
ably with an all-metallic design of sufficient mass that 
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RADIATION 


heat transfer by conduction produces the desired re- 
duction in stagnation area temperature. 


Concluding Remarks 


The calculations presented provide a ready means 
for estimating the structural temperature reductions 
attainable in certain radiation cooled areas of hyper- 
sonic glide vehicles operating at positive angles of 
attack. Temperature reductions are shown to be 
maximized when areas of intense convective heating 
are suitably insulated and when adjacent areas of lower 
heating are used as radiators of internally transferred 
heat. 

The structural temperature reductions are of suffi- 
cient magnitude to permit structural design of glide 
vehicles, operating at velocities up to and including 
orbital velocity, to be based on available high-temper- 
ature materials. Under certain design conditions, 
these temperature reductions may be accompanied by 
a reduction in total structural weight. 


Appendix 


In the analyses that follow, it is assumed that steady- 
state heating situations exist. In order to simplify the 
analyses, it is further assumed that the temperature, 
aerodynamic heating, and emissivity are constant over 


a given surface. 


Uninsulated Structure 

Case I. Wing Surface 
surfaces are treated as infinite gray-body areas with 
different emissivities and if internal reflection is per- 
mitted, the heat balance equations for the upper and 
lower surfaces, respectively, are 


If the wing upper and lower 


gu + oe*(T,! — Ty) = ceyTy' (A-1) 
dx = oe,T,' + oe*(T,* — Ty*) 

where 
e* = 1/|(1/evy) + (1/ez) — 1] (A-2) 


The temperatures of the upper and lower surfaces 
can be obtained from Eqs. (A-1) in the following dimen- 


sionless form 


Ty ade (e*/ev) + (ex/ev) [1 + (e*/ex) (Qu/qu) 4 
i. it + (e*/e,)|[1 + (e*/ey)] — (e*? ever) ) 
(A-3a) 
/ . « [1 + (e*/ev)] + (€*/ev)(qu/qz) ye 
Teo lo + (e*/ez)|[1 + (e*/ev)] — (e*? ever) S 
(A-3b) 


where by definition 
Tn = (9r/oer)"* 


If the emissivities of all surfaces are equal (i.e., 
éy = €, = ¢), then e* = €/(2 — e), and Egs. (A-3) re- 


duce to Eqs. (1) given in the text. 
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For the case of perfect internal heat transfer, both 
surfaces assume the same temperature which is deter- 
mined by the following heat balance equation 


gu + Qt = 2ocT* (A-4) 
or, in the dimensionless form, 
T/T.,1 = {(1/2)11 + (v/q) P44 (A-5) 


Case II. Leading Edge—For this case, the emissivi- 
ties of all surfaces are considered equal and sufficiently 
close to unity to neglect internal reflection. In order 
to provide a simple treatment of radiation interchange 
factors, the leading edge is considered infinitely long in 
the spanwise direction, which allows the assumption of 
plane radiation.° 

The heat balance equations for the areas involved 
are as follows: 
Stagnation area 


qsd + oelFy sly + ocd Fy gi w* + oelF, gl t' = 


(1 + Fs_y + Fs_w + Fs_1t)oedT 5‘ (A-6) 
Upper surface 
qul + oedF g yl s* + coed Fy yl w* + oelF, vl 1* = 
(1 + Fy Ss + Fy Ww + F, t)oedT y' (A-7) 
Closing web 
ocd F wis‘ + oelFy wiv — oelF, wi +* = 
(Fw S + Fy U + Fw toed T w'* (A-8) 
Lower surface 
ql + ced Fg ti 3° + oedF y L! w* + oelFy Li y* = 
(1 + Fr_es + Fr-w + Fr-v)oeT1* (A-9) 


Inasmuch as the four surfaces completely surround the 


enclosed volume, 


Fs ot Fs wt Fs t= 1 


Feat Fow+ Fou =10 (aig 
Fw-s + Fw-v + Fw-1. = 1 
Fi_s + Fr_w t+ Fr-v = 1 
From symmetry and the reciprocity theorem, 
Fy = Fi S» Fs U = Fs Ly Fy s = Fs a (A-11) 
Fy_t = Fr-v, Fw-v = Fw-r 
and dFs gg = lFy S» dFy gy = lFy Ww (A-12) 





48 JOURNAL OF THE AERO 

If Eqs. (A-10), (A-11), and (A-12) are substituted 
into Eqs. (A-6), (A-7), (A-8), and (A-9), the following 
set of equations, governing the temperature distribu- 


tion, is obtained: 

i 7.9 + 
Fs wl'y*| 

oF (d l) Fs u(7's* + 
Tw') + Fo_1T1'!| 

L') + Fs_wl's' 


(d L)Fs y(T 5* + 
") + Fy_1Tv"| 


T° = (1 Za sg! 4 Fs oll 


ie = CHS, 
(A-13) 


Tyw' — Fs o(Ty4 + 7 
2)(T.. 24 + 


where by definition 


T., s' = ds/ce 
T., u' = qu/ce 
I...‘ = Yi’ oe 


If the above set of Eqs. (A-13) is solved for 7s and 


use is made of interchange factor algebra which gives 
Fs w= 1-—2Fs_v | 

} (A-14) 

Fy, = 1 — 2(d/l)Fs_v\ 


7's; can be expressed as the following dimensionless 
function of #s_y and the ratio d/1: 


(T's Pa: s)* = gl[Fs U» (d !)| + 


1 = gl Fs U> (d !))} (ge ++ dz) 2¢s| (A-15) 
where 
g[Fs_v, (d/l) | = 0.5}0.75 + (d/l) F — {1.5+ 
(d/l) |\(d/1) F; eee 10.375 4 
[1.5 + 05D Fe.» [1.5 + 1.75(d/) + 
0.5(d/1)? 24 [(a/l) + 0.5(d/l)?|Fs_y*} (A-16) 





However, the interchange factor, Fs_,, can be ob- 
tained from the following equation which is derived 


from fundamental relations given in reference 5. 


Fs_y = (l/2d)[(d/l) —- V1 + d/)2 +1] (A-17) 
Therefore, the dimensionless temperature 7's5/7.,, 5 can 
be expressed as a function of d// only, or 
(Ts/Te, s)* = f(d/l) - 

[1 — f(d/l)|[(@v + gr)/2¢s] (A-18) 


For the case of perfect internal heat transfer, all four 
areas are at the same temperature which can be deter- 
mined from the heat balance of the entire leading-edge 
section. If all surfaces have the same emissivity, the 
heat balance equation is 


gsd + (gu + gz)l = oe(d + 21)T! (A-19) 


From Eq. (A-19), the following dimensionless temper- 
ature equation is obtained 


(T/T-, s)* = 
{1 + 2(1/d)[(qu + 


qz)/2¢s]}/[1 + 2(l/d)] (A-20) 
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Insulated Structure 

Case III. 
is insulated and the 
different, the following heat balance equations are ob- 
tained for the case where gray-body reflection exists 


Wing Surface —If the wing lower surface 


emissivilies of the surfaces are 
between two infinite wing surfaces: 
di, = gégl 9’ + 
(k t)(To —_ T;) oe*(T,' 


gu + oe*(T,' — Ty’) 


(k/t)(To — T;,z) 


a (- . (A-21) 


oe; T, . 


Eqs. (A-21) can be conveniently converted to the 
following dimensionless set which defines the temper- 


ature distribution: 


Ty 
‘iy oO 
|< ey) (Ty, i. o)' - (€o €; ACh te) i 
1+ (e* 
Ls oe ex = (A-22) 
=f -(- GO 
ee oO ( €0 Pe oO Pe oO aL 
t qh PZ (€0 e” ) ((To/7 e, o) o- (TT,  e o) | 
k T., 0 (TL T., o) jt — (Ty i a o)* 
where by definition 
T..o0 = (Gi Te) '/4 
If the emissivities are equal, ¢9 = €, oe = “<< 


e/(2 — e), and Eqs. (A-22) reduce to 


Ty = | = T..0)* = i (2 — €) (du al ; 
fae oO ins : aoa 


ti 4... 2 0% ) ty \} ! 
= 41i- _ A-23 
a: oO { a ie O Dx Oo f ' ?) 


(To/T.. o)} — (T2/T.. o)) 


2 — é€) 
\ . (TL tee)” — (Ty en 


kT, 0 

If all the emissivities are unity, or if the emissivities 
are large enough to assume that reflected radiation can 
be neglected, Eqs. (A-23) further reduce to Eqs. (7) 
given in the text. 
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Compressible Flat-Plate Boundary-Layer 
Flow With an Applied Magnetic Field’ 


WILLIAM B. BUSH* 
Space Technology Laboratories, Inc. 


Summary 


The laminar boundary-layer equations are formulated and 
solved for a flat plate in high-speed compressible air flow where 
equilibrium dissociation and ionization are assumed and where 
there is an applied magnetic field having its component normal 
to the plate proportional to 1/+/x. The skin-friction and heat- 
transfer characteristics are determined for free-stream velocities 
of up to 17,500 meters/sec. and magnetic fields of up to about 
1 weber/(meter)?. 

The results show that the skin friction and heat transfer at a 
given free-stream velocity decrease with increasing magnetic 
field strength, and the percentage reduction is constant along 
the length of the plate. They also exhibit the same hysteresis 
behavior as was first found in the case of magnetoaerodynamic 
Couette flow; however, for the flat plate the hysteresis effect 
disappears at a higher Mach Number. Furthermore, it was 
found that the reduction in heat transfer with increasing field 
strength is opposite in behavior from that for Couette flow. 


Symbols 

B = magnetic field strength 

B,, B, = Cartesian components of B 

Cy = specific heat at constant pressure 

E = electric field intensity 

Fy isc = viscous force 

G = shearing stress parameter 

h = specific enthalpy 

j = current density 

k = thermal conductivity 

M = Mach Number 

p = pressure 

P? = Prandtl Number (Pr = ne,/k) 

q = heat flux 

q = velocity 

Q = magnetic parameter 

Ri = viscous Reynolds Number 

t = independent variable 

i = absolute temperature 

ul = x and y components of velocity 

we = stream function parameter 

x, = right-hand Cartesian coordinate system 

z = space variable (s = v/v /x) 

a nondimensional product of density and viscosity 
(a = pn/pane) 

8 = nondimensional product of electrical conductivity 
and viscosity (8 = on/o7rnp) 

6 = boundary-layer thickness (or, in reference 3, dis- 
tance between walls ) 

c = heat-transfer parameter 

n = absolute viscosity 
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Me = magnetic permeability 

p = density 

o = electrical conductivity 

o7 = reference electrical conductivity (a7 = 100 mho/m.) 

T = local shearing stress 

PD\ ise = viscous dissipation 
Superscripts 

(> = nondimensional variable [see Section (I1)] 

( = nondimensional variable {see Section (IV )] 
Subscripts 

( Dr = values at base temperature 400°R. (= 222°K 

{ = reference value 

( ) = value at the plate 

( = value at the free stream 

( ) = quantity outside the boundary layer 

( )e = quantity inside the boundary layer 


(I) Introduction 


K )R HIGH VELOCITY FLOW of a gas past a body, the 
viscous heating in the boundary layer converts a 
large portion of kinetic energy into thermal energy. 
For extremely high velocities, this thermal energy is 
sufficient to partially dissociate the gas and even pro- 
duce a small, but often not negligible, degree of ioni- 
zation, and, as a consequence, the gas becomes an elec- 
trical conductor. Recently, studies have been made 
to see if the interaction of magnetic fields with these 
conducting gases will appreciably alter the skin fric- 
tion, heat transfer, etc., for such flows. 

One of the first such studies was made by Rossow, ! 
who obtained an approximate solution to the flat plate 
problem considering incompressible, laminar flow with 
constant electrical conductivity and a uniform mag- 
netic field applied normal to the plate. Lykoudis® 
went on to show that general ‘“‘similarity”’ solutions exist 
for the magnetoaerodynamic boundary layer for wedge- 
type flows if (1) the electrical conductivity (o) is con- 
sidered constant, (2) the induced magnetic field is as- 
sumed negligible in comparison with the applied mag- 
netic field acting perpendicular to the direction of flow, 
and (3) the intensity of this applied field varies accord- 
ing to a power law with respect to the distance along 
e.g., for a flat plate it is proportional to the 


the body 
Both of these papers show 


minus one-half power. 
that with increasing magnetic field both the skin fric- 
tion and heat transfer decrease. 

Bleviss,* appreciating the difficulties in including in 
the boundary-layer analysis the variation of thermo- 
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dynamic and transport properties, especially o, has 
investigated the magnetoaerodynamic hypersonic Cou- 
ette flow problem, since it contains many of the im- 
portant features of boundary-layer flow but at the same 
time can be solved exactly for all Mach and Reynolds 
Numbers with a minimum number of assumptions 
about the nature of the gas itself. His results for the 
nonzero heat-transfer case show that, for very high 
Mach Numbers, where the enthalpies at most points 
of the flow field are high enough so that ¢ varies slowly 
with enthalpy (4), the skin friction decreases mono- 
tonically with increasing magnetic field. For lower 
Mach Numbers, where the enthalpies at most points 
in the flow field are low enough so that o¢ varies rapidly 
with h, the skin friction as a function of magnetic field 
exhibits a hysteresis-like behavior. The skin friction 
decreases slightly as the magnetic field is built up, until 
a critical value of the field (B,,.) is reached, and then 
drops sharply to a much lower value. However, if 
the magnetic field is now decreased, the skin friction 
slowly increases until a second critical value of the field 
(B.,,) is reached (where B,,, < B,,,), and then jumps 
sharply back to a higher value. Bleviss predicts that 
this hysteresis effect will also be present in the flat- 
plate boundary layer. 

In this paper, an attempt has been made to see if 
Bleviss’ prediction can be verified. The problem 
considered is the flow of air over a flat plate (aligned 
in the positive x direction with x = y = O the leading 
edge). The temperature of the air outside the bound- 
ary layer is to be considered low enough so that the 
electrical conductivity in this region is negligible. The 
temperature of the plate itself is constant for all values 
For all points in the flow field, the air is con- 
Acting on this flow is a 


Ol %. 
sidered to be in equilibrium. 
magnetic field which (for mathematical convenience) 
has a y component proportional to x~'/?. This is the 
form of the magnetic field suggested by Lykoudis, but 
it was arrived at independently by the author. 


(IT) 


Equations of Motion 


For a nonferromagnetic substance such as air, the 
steady-state electrodynamic equations are 


curl B = yj (II-1) 
curl E = 0 (II-2) 
div B = 0 (II-3) 
j = o(E + qxB) (II-4) 
The mass, momentum, and energy equations are 
div (pq) = 0 (II-5) 
(1/2)p grad g? — pq x curl q + grad p = 
(jx B) + Fyisc. (II-6) 
pq:grad h — q-grad p = div (k grad 7) + 
P?/o + Pyise, (II-7) 


If the analysis is restricted to cases for which the 
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electromagnetic boundary conditions are such that E 
is zero, and the following nondimensional variables are 


introduced: 
“u = uv, B. = BoB. | 
o = 7'y, B, = B,*B, 
= pg" p, S = 27h, 
pb = p*p,V;? y = y*L, (11-8) 
ho h*V,? iv = u*u, 
T = T*V,?/€p, ¢ = o*¢, 
7 = n*n, k k*k, 


then for two-dimensional flow expressed in Cartesian 
coordinates, Eqs. (II-1)—(II-7) yield: 


(0/Ox*)(p*u*) + (0/Oy*)(p*v*) = 0 


p*u*(Ou*/Ox*) + p*v*(Ou*/Oy*) + (Op*/Ox*) = 
(n-/ pr VL) (0/Ox*)} n* [(4/3)(Ou*/Ox*) — 
(2/3)(Ov*/Oy*)]{ + (n-/p,V,L,) X 

(0/Oy*)} n*[(Ou*/Oy*) + (Ov*/Ox*) |} — 
(c,B,2L,/p,V,)o*B,*[u*B,* — 0*B,*] 


(II-9a) 


(II-9b) 


p*u*(Ov*/Ox*) + p*v*(Ov*/Oy*) + (Op*/dy*) 
(n,/ prV-L,)(0/Oy*)\ n* [(4/3)(Ov*/Oy*) — 
(2/3)(Ou*/Ox*) |} +(n-/p-V,L,) X 

(0/Ox*)} n*[(Ou*/Oy*) + (Ov*/Ox*) |} + 


(o,B,°L,/p,V,)o*B,*(u*B,* — v*B,*| (1-9) 


p*u*(Oh*/Ox*) + p*v*(Oh*/Oy*) — u*(Op*/dx*) — 
v*(Op*/Oy*) = (,/p,VrLr)(Rr/Cpnr) X 
) (0/Ox*) [R*(OT* /Ox*) | + 
(0/Oy*) [k*(OT*/Oy*) |} + (¢,B2L,/p,V)o* X 
[u*B,* — v*B,*|? + (n,/p,VL,)n*) — (2/3) X 
[(Ou*/Ox*) + (Ov*/Oy*) |? + 2(0u*/Ox*)? + 


2(dv* /Oy*)? + [(Ov*/Ox*) + (Ou*/Oy*) |?} (II-9d) 
(0B,*/Ox*) + (0B,*/dy*) = 0 (11-9e) 
(OB,*/dx*) — (OB,*/dy*) = 
(o,u,V,L,)o*u*(u*B,* — v*B,*| (I1-9f) 
If Le = (ogy V,) 
then nr/ Pr Veber = NO rhtr/ Pr = (Ro,) (11-10a) 
oB,?L,/p,V, = B,?/u,p,V,2 = Q, (II1-10b) 


and the other parameter, (k,/Cp,n,), is the reciprocal of 
the Prandtl Number. 

Since for air the logical reference magnetic per- 
meability is that for a vacuum, unless the temperature 
is extremely high or the density extremely low, it fol- 
lows that Rv, > 1. Further, if the flow is to be appre- 
ciably changed from what it would be in the case of no 
magnetic field, the magnetic forces must be of a magni- 
tude comparable to the primary fluid dynamic force. 
The magnetoaerodynamically important case, there- 
fore, is the one in which Q, is of order one. For tem- 
peratures up to at least 10,000°K., Prandtl Number is 
of O(1). 


(III) Boundary-Layer Approximation 


Since (Rv,)—! is small, one is led to a ‘‘boundary layer”’ 


or singular perturbation problem. Outside a_ thin 
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boundary layer adjacent to the body, in which the 
functions undergo rapid changes, the terms in the 
momentum and energy equations that are proportional 
to (Rv,)~! may be neglected. If, in addition, we fur- 
ther restrict the problem to one for which the temper- 
ature and density in the free streams are such that 
o* <1, the equations in the free stream become those 


for an inviscid, nonconducting fluid—namely, 


(O Ox1*) (pi *u1*) + (O Oy: *)(p:*21*) = () (III-la) 

pi *uy*(Ou,* 0x1") + pi *v,*(Ou,* Ov;*) + 
(0p,*/0x,*) = 0 (III-1b) 

pi*uy*(Ov,* 0x,*) + pi *2*(Ov,* Oy,;*) + 
(Op,*/Ov,*) = 0 (IIT-Ic) 

pi *u,*(Oh,* Ox,*) + pi*v*(Oh,* Oy, *) = 
U,*(Op;*/Oxy*) + %*(Op,*/Ov,*) (ITT-1d) 
(OB,,* Ox1*) + (OB,,* Ov;*) = 0 (III-le) 
(OB,,* Ox,*) a (OB,,* Oy,*) = 0 (III-1f) 


where the subscript | refers to the region outside the 
boundary layer. 

If the quantity 6 indicates the order of magnitude of 
the boundary-layer thickness (nondimensionalized in 
the same manner as x* and y*), the nonmagnetic 
boundary-layer arguments state that, in the layer, 


yo* = 0(6), u.* = O(1), 


v.* = 0(4), Rv, = 0(1/6*) 


co" = Ht). 
(III-2) 


where the subscript 2 refers to quantities in the bound- 
ary layer. 

Investigating Eq. (II-9e) for the boundary layer, 
since the applied field is one whose y component will 
be maximized, the nondimensional y component should 
be of order one, and, thus, the x» component would be 
expected to be at most of order (1/6). Substitution 
of these orders of magnitude into Eq. (II-9f) shows 
that (0B,*/Oy*) is of order one. By integration, it 
follows that, inside the boundary layer, B,* = B,*(x*). 
A further integration leads to the conclusion that B,* = 
B,*(x*) inside the boundary layer. In other words, 
the components of the field in the boundary layer may 
be approximated by their values at the plate. 

We shall restrict the analysis so as to consider the 
case where 


B,,* = K,*(x.*), Baz* = K.*(x.*) (IIl-3) 


The nondimensionalized boundary-layer equations, 
taking into account the preceding orders of magni- 
tude, thus become 


(O Ox2*) ( po*uto*) + (Oo Ovyo*) ( p2*v*) = (0 (III-4a) 


po*uo* (Ous* /Ox2*) + po*v2*(Ou2*/Ov2*) + (dpo*/dx2*) = 
(Rv,) —1(0/Oy2*) [q2* (Ou2* /Oy2*) | — Q,02*u2*(B,2*)? 
(III-4b) 


Op.*/Oy.* = O(1) (III-4c) 
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po*uo* (Oh»* Ox2*) + p»*v2*(Oh.* Oye*) = 
u2*(dp.*/dx.*) = (Rv,Pr,)—' X 
ko* (OT >* Ovye*) | + (Rv,) lno* x 
(Om2* /Oy2*)? + O,02*(u2*B,*)? 


( Oo Ov.* ) 





(III-4d) 


Derivation of the Equations Suitable for 
Machine Computation 


(IV) 


From the previous section it can be seen that the 


dimensional laminar boundary-layer equations for 
compressible flow over a flat plate including magnetic 


terms are 


(0/Ox)(pu) + (0/Oy)(pr) 0 (IV-la) 

pu(Ou/Ox) + pv(Ou/Oy) 
(0/dy)[n(Ou/dy)] — ouK,2 (IV-1b) 

pu(Oh/Ox) + pv(Oh/Oy) 
(0/Ov) [(n/Pr)(Oh/Ov) | + n(Ou/dyv)? + ou?K,? (IV-Ic) 


‘In order to reduce these partial differential equa- 
tions to ordinary differential equations, the nonmag- 
netic analysis of Hantzsche and Wendt‘ and Moore® 


is used. In this analysis it is assumed that 
u u(z), Ah h(z) (IV-2) 
where z (y/~/x). Since, in general, the density 


and transport properties can be expressed as functions 
of enthalpy and pressure, and, for flat plate flow, the 
pressure is a constant, it follows that 


p = p(z), 9 = (2), Pr Pr(s), o 


For the magnetic terms to be incorporated into this 
analysis examination of the equations shows that A, 
must be 


K, = (BoVLy)/V'x (IV-4) 
where (By+/L,) is a constant. 

This magnetic field configuration has a singularity at 
x = 0, and such a field would not be physically realiz- 
able. 
layer equations for a nonconducting fluid, we know that 


However, from previous studies of the boundary- 


the solutions cannot be expected to be a good approxima- 
tion near the leading edge—e.g., for the nonmagnetic 
solution of Hantzsche and Wendt, the normal com- 
ponent of the velocity, v, is proportional to x~'/? 
Therefore, we see that the assumption as to the con- 
figuration of the magnetic field breaks down in the 
same region as do the boundary-layer equations, and 
is valid in the region where the boundary-layer equa 
tions are valid. 

Transforming the independent variable from z to the 
velocity variable u and introducing a function G(u) = 


n(du/dz), the combination of the momentum and 
energy equations yields 
(dh/du)(dG/du) = (d/du)|(G/Pr)(dh/du) | + 

G + [onu(Boy*Lo)/G\[(dh/du) + u}] (1V-5) 


Combining the mass and momentum equations gives 
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G(d°G/du*) = —(pnu/2) + 


G(Bo?Lo) (d/du)(anu/G) (1V-6) 


In order that the analysis of Moore may be followed 
as closely as possible, his method of nondimensionaliz- 
ing the variables is used: 


Velocity: fi = u/Vhp (1V-7a) 
Shearing Stress: 

G(a) = G(u) V (1 2) ppnp(hp)*!? (IV-7b) 
Enthalpy: h = h/hp (IV-7c) 


where the subscript B indicates the value of the air 

properties at the reference temperature of 400°R. 
(= 222°K.). Eqs. (IV-5) and (IV-6) then become 

(dh/dit)(dG/di) ill (d dit) [(G, Pr)(dh/di)| + 

G + Q(8u/G)((dh/da) + a] 

G(d°G/di?) = —ati + OG(d/di)(Bi/G) (IV-9) 


(IV-S) 


where 


a = (pn/prnn) = a(h) (1V-10a) 


8 = (on/ornn) = B(h) (IV-10b) 


O = (207rBy*Lo pe hs) = constant (IV-10c) 


and or is the reference electrical conductivity and is 
taken to be 100 mho/m. 

Thus, Pr, a, and 8 contain all the physical properties 
of air including the effects of dissociation and ionization. 
These quantities, as mentioned previously, can be ex- 


TABLE 1 
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pressed as functions of / or h, which is also affected by 
dissociation and ionization. 

Moore makes one more transformation to eliminate 
the infinity value at the upper limit (1/G —~ ©) which 
is not suitable for machine computation. This trans- 
formation is 


G = di/dt (IV-11) 


A new variable w, is also introduced and defined by 


dw/dt = at, w(0) = 0 (IV-12) 
The two equations then become 
dG/dt = —wG + Opa (IV-13) 
@h/dt? = —wPr (dh/dt) + ' 7 
(1/Pr)(d Pr/dh)(dh/dt)? — PrG? — PrQpi*  (IV-14) 


Thus, the equations in the form required for the digital 
computer are 


di/dt = G (IV-15a) 


dw/dt = att (1V-15b) 
dG/dt = —wG + Qpi (IV-15c) 
dh/dt =¢ (IV-15d) 
d¢/dt = —wPr¢ + (1/Pr)(d Pr/dhjge? — 


PrG? — PrOBa? (IV-15e) 


where the original independent variables, x and y, are 
related to ¢ by 


*/ 
y = (2nex/ppV hp)'” | (n/ne)dt (1V-16) 


so that when y = 0, ¢ = 0. 
The boundary conditions at the wall for these equa- 


tions are 


u(O) = O (IV-17a) 
w(0) = 0 (IV-17b) 
G(0) = G, (IV-17e) 
h(O) = hy, a specified value (1V-17d) 
(0) = Su (IV-17e) 


The values G, and ¢, must be adjusted so that as ¢ 
approaches “‘infinity’’ # and / approach their free- 


stream values (and (@ approaches zero). 


The Variation of h, a, and B with 7 for a Pressure of 1 Atm. 7 PABLE 2 
eee poeennonaiases — The Variation of 2, a, and 8 with 7 for a Pressure of 10~* Atm. 
si ag a h a 8 LN ee oe Satadecree sil 
= ‘ 6 SALT mn” s..) h a B 
200 0.88 1.03 * 10° = i a a 
300 1.35 9.46 X 107! 200 0.88 1.03 X 106° 
400 1.79 8.81 < 107 300 1.35 9.46 KX 107 
500 2.22 $.21 «% 10- 100 1.79 8.81 <X 107 
1,000 4.65 6.39 * 107! 500 2.23 $.21 * 1 
1, 500 7.29 5.39 & 107 1,000 4.65 6.42 X 10™ 
2? 000 10.2 5.00 & 107 2,000 10.6 5.05 X 107 
3,000 18.2 4.22 = 107 '.36 x 107 3,000 31.4 3.60 KX 107 1.74 X 10 
+, 000 34.4 8.32 107 1.74 X 107! 4,000 42.4 2.82 * 107 9.10 K 107! 
5,000 45.9 2.90 * 107} 2.14 & 10° 5,000 106.0 1.98 X 107 6.14 X 10° 
6,000 67.0 2.56 < 1073 1.07 * 10! 6,000 168.0 1.42 <X 107 3.64 & 10! 
7,000 117.0 2.14 x 107 3.54 X 10! 7,000 185.0 lad & 10™ 1.22 X 102 
8,000 170.0 1.89 X 107 1.14 * 10? 8,000 221.0 1.42 x 10 2.07 & 102 
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(V) The Thermodynamic and Transport 
Properties of Air Used for the Calculations 


As has been seen in Section (IV), to solve the flat- 
plate boundary-layer equations with a magnetic field 
applied, the variables Pr, a = (pn/psn,), and B = 
(on/orns) must be found as functions of the dimen- 
sionless enthalpy, 4. For the present work they were 
found for pressures from 1 to 10~* atmospheres and for 
a temperature range of 200° to 8,000°K. for air in 
thermodynamic equilibrium. Constant pressure en- 
thalpy and density variation with temperature were 
obtained from Gilmore’s report,® while the temper- 
ature variation of viscosity was derived by Green’ and 
Greifinger.2 The behavior of electrical conductivity 
with temperature is based on the work of Lin and 
Lamb,’ Meyer," Bleviss, and the author.'! Reference 
8 indicates that, for the temperature range of interest, 

Pr = 0.7 = constant (V-1) 
is accurate enough for the purposes of this paper. 

a, B, and h as functions of temperature for pressures 
of 1 and 10~* atm. are tabulated in Tables 1 and 2. 
Fig. 1 shows a as a function of h for pressures of 1 and 
10-* atm. Fig. 2 shows @ as a function of / for pres- 
sures from 1 to 10~* atm. 

Due to the difficulty in using graphical data for solv- 
ing equations on the digital computer, it was thought 
to be advisable to approximate the behavior of a(/) 
and @(h) by functions more easily absorbed into the 
digital computer’s routine. 

Fig. 1 shows that a(/) is almost totally insensitive 
to wide variations in pressures and can be approximated 
quite well by 


(V-2) 


a= jr —0.35 
This approximation is also shown in Fig. 1. It is seen 
that this is a better fit for 10~* atm. than for 1 atm. of 
pressure. 

Fig. 2 shows that @(h) has a definite but relatively 
weak variation with pressure so that if the problem 
were to be carried out exactly, a solution would have to 
be found for each pressure. However, from the data 
it can be seen that a good approximation is of the form 

B = Boh" (V-3) 
where the pressure dependence shows up almost en- 
tirely in determining the constant, 8, and has very 
little effect in determining the exponent, x Further- 
more, since 8 never appears alone in the problem but 
always appears in the product (Q@), the variation of 8 
with pressure can be absorbed in the parameter Q. 
For this paper 8 was taken to be 


8 = 1.20 X 10-%*-7° for h > 4.65 (T > 1,000°K.)\ 
= 0 for h < 4.65 f 
(V-4) 
This approximation is shown in Fig. 2. It can be seen 
that this approximation fits the results for p = 107% 
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Electrical conductivity-viscosity parameter vs. enthalpy 
parameter 


Fic. 2. 


atm. best, but, with proper evaluation of Q, can serve 
for all pressures. 
Specifically, for p = 1 atm., 

By ~ 5.68 X 10-3 (V-5) 
that is, at 1 atm. 6» is 47.4 times as large as 8, for 10 
atm., and Q at | atm. needs to be only (1 47.4) the size 
of Q at 10~-* atm. to produce the same magnetic effect. 
Therefore, to convert the results found for p = 10 
atm. to those for 1 atm., divide Q by 47.4. 


Since 
O = 207B,"*L osV hea = 2o7B,*x/ ppVhp (V-6) 
it follows that, for (%)14 = (X),0-a4, 
(By)1a/(By)10-14 = ((Opp) 14 (O Pp)10-sa |!” 
10920 (O)14/(O)i-s4 (V-7) 
and to produce the same magnetic effect 
(B,):a/(B,)w-sa * V108/47.4 = 4.61 (V-8) 


Discussion of the Computing-Machine 
Solutions 


(VI) 


The five first-order differential equations, Eq. (IV- 
15), modified by Eqs. (V-1) and (V-4) (1.e., 6 = 10 
atm.), were solved by assigning the boundary condi- 
tions at the plate (¢ = 0), Eq. (IV-17), and allowing 
the machine to proceed with the solution as ¢ goes to 
infinity. The values of hf; and a, the values at the 
outer edge of the boundary layer, were thus deter- 
mined from the solution. For all cases, /,, was taken 
to be 10, which corresponds approximately with a tem- 
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perature at the flat plate of 2,000°K., as higher wall 
temperatures were felt to be unrealistic. 

Specifically, a series of runs was made in the follow- 
ing manner: a value of G,, was selected and held fixed 
while ¢,, was allowed to vary over a large range and the 
resulting values of 7; and h; found; when the range of 
¢» was sufficient to yield the required information, a 
new value of G, was selected and the procedure was 
repeated. This scheme was carried out for a wide 
range of the parameter Q. Unfortunately, due to the 
large number of runs required, to reduce machine time, 


oO 
= 





Fic. 3. Shearing stress parameter at the plate vs. velocity 
parameter outside the boundary layer forh; = 1, hw = 10,p = 


10-2 atm. 
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p = 10°% atm 


the numerical method used did not compute the bound- 
ary-layer profiles of velocity and enthalpy. 

When the runs were completed, for each value of Q 
there were two graphs, one graph showing /; vs. ¢,, for 
varying G,. and the other showing a; vs. ¢, for varying 
G,. Cuts were taken in the first graph for a given h; so 
that G,, and ¢» were determined for this given /;. 
These values of G,, and ¢» then determined 7%; from the 
second graph. 

In Figs. 3, 4, and 5, G.. and ¢,, are plotted as functions 
of a; for hs = 1 and O = 0, 1, 10, and 100 (with /,, = 
10 and » = 10~ atm.). 
functions with i; = 10. 


Figs. 6, 7, and 8 show these 


Since 


tT» = [n(0u/dy) |. = [pane(he)*!?2/2x |G, (VI-1) 


Gow = [R(OT/Oy) |o = Prlpane(hp)*!?/2x]'/"¢, (VI-2) 


it follows that, for a given x, 7, and p at the wall, 


(Te)1 (Tw)2 = (Gy) (Gin)2 (VI-3) 


(Gw)1 (Gwe _ ey? (Cw) (VI-4) 


Using these last two equations, it was possible to use the 
preceding data to find [t./(t)war] and [Gw/(Gw)na 
functions of Q for given M,; (or iis), hy, hs, and p. Fig. 
9 shows these relationships for 1J; = 25, and h; = 1 
(and h, = 10, pb = 10-* atm.). Fig. 
for M; = 30 and h; = 1. 
h; = 10, but in the former /; = 7.5, while in the latter 
M; = 10. 

From Figs. 3 through 12 an important fact stands 
out, that the skin-friction and heat-transfer curves are 
divided into three distinct regions by the vertical tan- 


as 








10 shows them 
Fig. 11 and Fig. 12 are for 


gents to the curves. In region 1, to the left of the 
tangents, the flow in the free stream has such a low 
velocity that the temperatures in the boundary layer, 
regardless of the value of Q, never are high enough to 
produce significant ionization; hence, the fluid behaves 
like a nonconductor and there is no magnetic effect. In 
region 3, to the right of the tangents, the free-stream 
velocity is high enough so that there is very significant 
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ionization in the boundary layer, and in this region the 
magnetic field has a large effect and reduces the skin 
friction and heat transfer to a very small fraction of their 
corresponding nonmagnetic values. 

In region 2, between the tangents, there are three 
values of G, and ¢,, for every value of #;. It is clear 
that this middle region represents an unstable flow 
condition. This means that if the velocity is increased 
from zero, for a given value of QO, G,, (or ¢,) will follow 
the upper stable branch until the right vertical tangent 
is reached. If the velocity is further increased, G,, (or 
¢») will drop abruptly to the corresponding value of 
tis on the lower stable branch. If, however, the ve- 
locity is now decreased, G,, (or €») will continue along 
the lower stable curve until the left vertical tangent is 
reached, and, for a further decrease in velocity, will 
jump back to the upper stable (nonmagnetic) branch. 
The boundaries for these regions are shown in Figs. 13 
and 14. Thus, for O < 100 (which is just about the 
practical upper limit on the magnetic parameter), if 
M; is less than 21.4 and 6.4 for hs; of 1 and 10, respec- 
tively, the skin friction and heat transfer are unaffected 
by the application of a magnetic field. 

At this point, it is appropriate to see just what mag- 
netic field strength must be applied to produce the 
different values of 0. From Section (IV), Q is defined 
as 


Q = 2a7rBy7Lo pa kp (IV-10c) 


where o7 = 100 mho/m. Since, for a pressure of 10~* 
atm., pp = 1.59 X 107* kg./m.* and V he = 472 m. 
sec., it follows that 


O/By?Ly = 2.66 X 102[(m.-sec.)/(ohm-kg.)] (VI-5) 


Thus, for Q = 100, remembering that B, at the wall is 
equal to K, = Bo(Lo/x)?, 


BLy = B,?-x = 0.376[(ohm-kg.)/(m.-sec.)] (VI-6) 
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For x = 1 meter, this becomes 


B, = 0.615(weber/m.*) = 6.15 & 10% gauss (VI-7) 


For the same pressure and x, it also follows that for 


0 =1, 


B, = 6.15 X 10? gauss (VI-S) 


For the same x and Q, but p = 1 atm., B, must be in- 
creased by 10*/?. 

Because of the differences in the formulation of the 
Couette flow and boundary-layer problems, it is quite 
difficult to make quantitative comparisons of the two 
However, it is profitable to make a quali 
The most im- 


theories. 
tative comparison of the two theories. 
portant facts found about the similarities and dissimi 
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larities of the boundary-layer theory and the Couette 
flow theory are as follows. 

(1) The boundary-layer theory predicts that the 
heat transfer for a given Mach Number decreases with 
increasing magnetic field strength, which in general, 
is just the opposite behavior from that found in the 
Couette flow theory. This difference in the heat 
transfer between boundary-layer theory and Couette 
flow is just what was expected by Bleviss. However, 
rather than the 15 to 30 per cent decrease expected 
by Bleviss, decreases of more than 80 per cent have 
been found for the field configuration of this paper. 

(2) In reference 3, the product of the normal mag- 
netic strength at the wall, B,, and the boundary-layer 
thickness, 6 (for reference 3 the distance between the 
moving and the stationary walls), rather than the 
quantity Q was the measure of magnetoaerodynamic 
effect on the flat-plate boundary-layer flow. For this 
reason it appears profitable to determine the product 
B,é as a function of Q for our problem. 

Since 


t 
= (2npx pnV hp)? ‘ff (n np)dt (IV-16) 
0 
then 
bt 
6 = (2nex onV hp)! ‘ff (n/np)dt (VI-9) 
and, thus, 


ét 
B = ByV x(2np ppv hp) ‘ff (n np)dt 
0 


a 


ot 
= (207By"Lo paV hp)"(ne or)’ ‘ff (n/np)dt 
0 


bt 
= Q"?(np/or)' “f (n/np)dt 
(VI-10) 


Since ng ~ 1.5 X 10° kg./m.-sec. and a7 = 100 mho/ 
m., the quantity (ns/or)/* ~ 4 X 10~* weber/m. = 4 
Further, 6, ~ 2 and the average value of 
5, so Eq. (VI-10) can 


gauss-m. 
(n/ng) in the boundary layer ~ 
be written as 


B,6 ~ 400"? gauss-meters (VI-11) 


From Figs. 9-12, it is seen that the magnetoaerody- 
namic effects are first really important at Q ~ 100, or 
for B,5 ~ 400 gauss-m. This is roughly the magni- 
tude at which the magnetoaerodynamic effects become 
important for Couette flow. 

(3) The hysteresis effect disappears at a much 
higher Mach Number for the boundary-layer theory. 

It is important to note that, in addition to neglecting 
the induced magnetic pressure, no mention is made of 
the finite mean free path (or molecular nature) of the 
fluid and, hence, the possibility of Larmor motion is 
ignored. This means that the conductivity is con- 
sidered unaffected by the magnetic field. This ap- 
proximation is good when (t,/t,) < 1 where ¢, is the 
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effective electron collision time and ft, is the Larmor 
period. 

For the magnetic field, B,, less than 0.5 (webers/m.’), 
the effective electron collision time is 


te ~ 10~15(po/ pp) sec. (VI-12) 
The Larmor period is 
t, ~ (10-"!/B,) sec. (VI-13) 
if B, is measured in (webers/m.’). Thus, 
(t./tr) ~ 10-*(p0/ps)B, (VI-14) 
Since 
B, = V 0/x(pp ‘po)'/2  (weber/m.?) (VI-15) 
then 
(te/tr.) ~ 10-°VQ/x(po/pz)!? — (VI-16) 


If it is required that (¢,/t,) < 0.1, it follows that 
O < 100 x (ps/po) (V1-17) 


For x = 1 m., 


Q < 100(px/ po) (VI-18) 


Therefore, for p = 1 atm., the conductivity should not 
be affected by the magnetic field if Q < 100. How- 
ever, for a pressure of 10~* atmosphere Q must be less 
than 0.1, if the conductivity is not to be affected. 

If the properties of the conductivity are not drasti- 
cally changed for (¢,/t,) ~ 1, then the present results 
should be applicable for 
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O = 10,000 x (px/ po) (VI-19) 
This means that for x = 1 m., we may consider 0 = 


10,000 as an upper limit for p = 1 atm. and 0 = 10 as 
an upper limit for p = 10~* atm. 

Finally, it should be stated that in this paper many 
simplifications have been included for mathematical 
convenience in an attempt to understand the magneto- 
aerodynamic phenomena in the flat-plate boundary 
layer. A more rigorous analysis, among other things, 
would either more fully justify the assumption of ther- 
modynamic equilibrium or discard it. Also required 
would be a more complete knowledge of the thermody- 
namic and transport properties at high temperatures. 
Thirdly, such an analysis would include taking into 
account different magnetic field configurations. 
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Temperature Profiles in a Finite Solid With 
Moving Boundary’ 


C. FORBES DEWEY, JR.,* STEWART I. SCHLESINGER,** anp 
LAWRENCE SASHKIN*** 


Aeronutronic, a Dwiston of Ford Motor Company 


Summary 


A numerical solution is presented to the transient heat con- 
duction equation for a cylinder of finite thickness with one 
moving boundary. The implicit method of solution is developed 
with conductivity as an arbitrary function of temperature. 
Application is made to a sample case of re-entry heating encoun- 
tered by aerodynamic bodies, with erosion by sublimation and 
combustion occurring at the body surface. 


Symbols 
a = original radius of moving boundary, ft. 
b = fixed boundary radius, ft. 
r = thermal capacitance B.T.U./lb.-°R. 
g, J = dimensional constants 
K = thermal conductivity, B.T.U./sec.-ft.-°R. 
N = number of space divisions in the slab 
q = heat flux at the boundary, B.T.U./ft.?-sec. 
r = radial coordinate, ft. 
Ss = radius of moving boundary at any time, ft. 
t = time coordinate, sec. 
T = temperature, °R. 
uy = free-stream velocity, ft./sec. 
v = transformed dependent variable 
z = nondimensional wall thickness 
p = material density, lb. /ft.* 


Subscripts and Superscripts 


b = reference to stationary wall 

z = space-point reference, 2 = iAz 

j = time-point reference, t = jA/ 

m = order of the approximation of the temperature profile 


reference to initial (¢ = 0) value 
reference to moving wall 


0 = 


AY = 


Introduction 


1822, 


WORK IN the 
in solid bodies 


Sn FOURIER’S CLASSICAL 
transient temperature distribution 
has been an important problem for 
mathematicians. Carslaw and Jaeger' present analyt- 
ical solutions to a wide variety of boundary value 
problems for both the infinite and semi-infinite body. 
Many other analytical solutions, graphical methods, 
and numerical solutions are available in the literature. 
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Because of the mathematical difficulties involved, 
much less work has been done considering variable 
conductivity.2:* For cases where thermal stress is 
important, or where physical properties depend on 
temperature, variable conductivity must be considered. 
Unfortunately, the effects introduced by variable con- 
ductivity are nonlinear, and a change in boundary con- 
ditions—e.g., heat flux or wall temperature, cannot be 
accounted for by using the usual methods of solution. 

Another complexity to the transient heat-conduc- 
tion problem which has been receiving an increased 
amount of attention is the case of a moving bound- 
ary.‘~* In particular, re-entry body surfaces subject 
to ablation or sublimation would fall within this cate- 
gory. Certain of these cases with a semi-infinite body 
may be solved analytically; but, in general, the time- 
dependence of the boundary location makes the result- 
ing differential equation difficult, if not impossible, to 
solve for arbitrary boundary heat fluxes, erosion rates, 
and initial temperature distributions. 

Through the use of high-speed computing machines, 
a general numerical solution to the combined problem 
of heat conduction in a finite solid with moving bound- 
ary and variable conductivity becomes feasible. In 
this paper, an implicit finite difference technique is used 
which permits arbitrary variation of the temperature 
or heat-flux boundary conditions with time. For use in 
problems such as aerodynamic heating where the heat 
conduction solution in the solid must be matched with 
external chemical kinetics, the finite difference scheme 
lends itself to an iterative matching of either heat flux or 
temperature. The initial conditions may also be made 
an arbitrary function of thickness. 


(1) The Heat-Conduction Equation 


(1.1) Transformation of the Differential Equation 

In this paper, consideration is given to the solution of 
the pure radial heat-conduction problem with variable 
conductivity which is governed by the differential 
equation 


(O7'/ot) = (1) 


(1/per) (0/0r)[Kr(OT Or) | 


The boundary conditions for radial conduction are 
given in terms of the heat flux or temperature at the 
two walls; these possible boundary conditions are: 


oT /or = —(1/K)q.(t), rv = s(t) 


oT /dor = —(1/K)q(t), r = 0 
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T=T,;(t), r = s(t) (4) 
T=T,(t), r=) (5) 


An arbitrary temperature profile is given as the initial 
condition 


f=, t= s =e 


The radius, s, of the moving boundary has been speci- 
fied as a function of time. This is done for computa- 
tional purposes only, since the motion of the boundary 
may easily be specified as a function of temperature or 
heat flux at the surface. In such a situation, the tem- 
perature (or heat flux) used to compute the boundary 
location is taken as the time average over the compu- 
tation interval, using the previous iterative step in the 
solution. 

To eliminate the moving boundary, a nondimensional 
space variable is introduced :* 


z= (r — s)/(b — s) (6) 


Because of the necessity to obtain a “‘smooth’”’ behavior 
of the final difference equations near the boundaries, 
we also make a transformation of the dependent vari- 


1 
7 =| K(r) dr (7) 
TR 


The transformed differential equation and boundary 
conditions may then be written: 


able: 


Ov /Ot = [K/pc(b — s)?](0?v/dz7) + 
({(1 — 2)/(b — s)](ds/dt) + 
}(K/pc)/[2(b — s)? + s(b — s) ]} )(de Oz) (8) 


(Ov/Oz) = —(b — s)q,(t), 2 = 0 (9) 

(Ov/Oz) = —(b — s)m(t), z= 1 (10) 

v=u,(t), 2=0 (11) 

v= (t), s=1 (12) 

The initial values are now written v(z) for ¢ = 0 and 


$ = @. 
Eqs. (8)—(12) completely define the problem under 
consideration. 


(1.2) Derivation of Finite Difference Relation 
Using time centered differences to replace all deriv- 
atives, the following relations are obtained: 


Ov/Ot == (V; 541 — 0; 5) /Al 

Ov Oz ~ (Vi4a, 541 — Dy 1, 42 es, 9 = Dsa4, $) 4 Az 
9 4 ‘ 

07/02? = (Vis, p41 — 204, 541 + V1, p41 + 


26 e ‘ g2 
Te a a 201. 5 +- Vi-1, j) 2Az 


where 7 represents the space index and j the time index. 

* An analogous solution to the one-dimensional slab with con- 
stant conductivity and moving boundary was solved by Ehrlich,” 
using the physical coordinate system rather than the transformed 
z-space defined here. For his case, Ehrlich defines the appro- 
priate method of dealing with boundary grid points. The trans- 
formed space variable, used by Otis,’ was also adopted by Murray 
and Landis" in treating the two-phase melting problem. 


JANUARY, 1960 


Substitution into Eq. (8) yields the difference equation 


Vi, j41 = U1, ; + [KAt/2pc(b — s)*Az] X 

Was st = 20;, j+l -+- Vi-1, j+1 + v4.5 = 
20;, 5 + vi-1, 3) + (At/4Az) ([(1 — 2)/(6 — s)] X 
(ds/dt) + }(K/pc)/[2(b — s)? + s(b — s)]}) X 


Wat, pt = tee pe Se 8a, 3 — Vi) «= (88) 


The above difference relation has been linearized in 7 
by assuming that the coefficients multiplying the inde- 
pendent variable are known functions of space and 
time. Insofar as they depend on the temperature, they 
are evaluated for the time interval using the temper- 
ature profile obtained from the previous step of an iter- 
ative process averaged with the known solution at the 
beginning of the interval. The temperature profile is 
iterated by successive substitution to the desired degree 
of convergence before continuing to the next time step. 

The difference relations for the temperatures at 


time ¢ = (j + 1)4¢ are written in the following form 
for ease of computation: 
Bovo + Cov1 - Dy / 

M-1+ Bo, + Cait1 = Di, 1<1< N —1> (14) 
Ayty-1 + Byty = Dy 


where the coefficients are given below. 
Given External Wall Temperature (i = 0) 

By = 1, CGo=0, Do = 2, j4:(0) (15) 
Given External Wall Heat Flux (4 = 0) 


Bo = 1 + 2B, Co — — 2Bo { 


; ‘16 
Dy = Vo (1 _— 2B) a 2B0Vi, j = Oo » 
where Qo = volvo — Be)gs(t) 
For all space points within the solid, 
A,= -Bi+y, Bi =1+28, CG: = —Bi - 7:| 
Di = 14,5 + Bivins, 5 — 201, 5 + Vi, —) + 
VilVini, 5 — U% wa 
(17) 
Given Inner Wall Temperature (1 = N) 
An = 0, By = L, Dy = Up i41(0) (18) 
Given Inner Wall Heat Flux (4 = N) 
Ay > —2Bn, By =1+ 2By 
Dy = 0p (1 — 2By) + 2Bntw-1, 53 — On (19) 
Qn = unlynw + By) @(t) \ 
where 
B; = [KAt/2pcAz(b — s)?] O<i<sN 
vi = (At/4Az) ([(1 — 2)/(b — s)](ds/dt) + 
}(K/pc)/[2(b — s)? +5(b — s]}}), O<i<WN 
a; = 442 (6b — ss), allz 


The heat-flux boundary conditions of Eqs. (16) and 
(19) were obtained by using the standard three-point 
approximation for (0v/0z) and eliminating the ficticious 
temperatures ati = —1l1 andi = N + 1 by combining 
Eqs. (9), (10), and (13). As pointed out in reference 
18, it may be preferable in some cases to include the 
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TEMPERATURE PROFILES IN A 


second term of the Taylor expansion at the boundaries. 
The transformation given in Eq. (7) reduces the possi- 
bility of large computational errors being introduced at 
the boundaries due to the variation of conductivity 
with temperature. 

A special method has been developed to handle the 
series of simultaneous linear relations which appear in 
Eq. (14).° This method is the equivalent of Gaussian 
elimination, where the constants 
eliminated in a “forward” direction and determined by 


are successively 
substitution in a “rearward’’ direction. 

For the relations of Eq. (14) which are linear in the 
undetermined quantities v;, we define 


W = B, W,=B,-Abis, 1<i<N 
b =C/W, OSi1<SN—-1, f= Do/We 
g = (D,—Agia)/W, 1<i<N 


(20) 


From the above definitions, the constants are succes- 


sively eliminated in the direction of increasing 7. The 
values of v; are found from the relations 
if — 3) 2K At . , &lAz 1 — sds 
= = -{ —sin? + 
(¢ + 1) pcAz?(b — s)? 2 b — sdt 
tlAs tlAs tlAz 
— “ . ‘ S$ ~ S$ ~ 
- —e —— +4 h sin cos 


where p = (2KAt)/[pcAz?(b — s)*] 


_ 
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-~_ ft, = f= b; Vi4+ty 0 < t < N- 1 (21) 
This elimination scheme facilitates the use of high-speed 
computing machinery by reducing the time and mem- 
ory capacity necessary for solution. 
(1.3) Stability Analysis 
Consider the difference equation 
({((1 — 2) /(b — s)](ds/dt) + 
 (K ‘pc) ‘[2(b — s)? + s(b — s)}}) y 
(Vi4a, j41 a 1, j+1 + Vi4i, 7 — 04-1, j) x 
(At 4Az) + [KAt/2pcAz*(b — s)?] X 
=, 20, +1 + U4 1, f+) TT Vesa, 3 _ 
))) 


24;, i + Ui, ;) (22 


OT == §) a 


Vi, j+1 Vi, = 


(Dini, j41 


from the point of view of stability. It may be shown 
that the truncation error €,, » satisfies the same differ- 
ence relation as v;,; for a system linear in »;, ;. Use 
of the usual Fourier decomposition of error employing 
Azonm! Vields the following equation for the 
e™. 


Com  « 


error growth ¢ = 


K / pc At/,. . &As ElAz 
- (7 sin cos 
z(b — s)? + s(b — s) |} Az 2 2 


(23) 


(b — s)|(ds/dt) + \(K/pc)/[2(6 — s)? + s(b — s)]})(At/ Az) 


h = ({(1 - 2) 

. 1 — psin?(t/Az/2) + 1h sin(élAz/2) cos (€/Az/2) 
— + p sin?(t/Az/2) — ih sin (é/Az/2) cos (Az /2) 
i [1 — p sin?(t/Az/2)]? + h? sin?(¢/Az/2) cos?(¢lAz/2) 

ic? = — a 


The error growth/time step, ¢, thus satisfies the rela- 


tion 
|< | < lforall p > Oandh (24) 

In this derivation of the stability criterion, it has been 
assumed that the conductivity and boundary moving 
rate are slowly varying functions of time. With these 
assumptions, it can be concluded that the difference 
equation is unconditionally stable. 

The selection of this particular method for the solu- 
tion of the partial differential equation was due to the 
great stability of this numerical integration procedure. 
Unlike the explicit difference procedures in which all 
values of temperature at the new time are described 
explicitly in terms of the temperature distribution at the 
preceding time,'* it is not necessary to integrate using 
small increments in the time variable to retain sta- 
bility—i.e., to prevent the exponential growth of errors 
from step to step. 

Using this implicit procedure, the integration interval 
is selected simply to control truncation error in each 
step; one need not be concerned with the magnification 


-of this error in successive steps of the integration. 


[1 + p sin? (/Az/2)]? + h? sin? (/Az/2) cos*(tlAz/2) 


Such implicit integration procedures as these have been 
used with tremendous success in such diverse applica- 
tions as reactor design calculations, hydrodynamics, 
and diffusion of fluids through porous media. 


(2) 

The difference equations presented above may be 
applied to re-entry bodies experiencing surface erosion 
due to combustion and ablation. In this application, 
the rate of removal of the surface material depends upon 
the wall temperature and the external flow conditions. 
In turn, the wall temperature is determined by the flow 
field (the net heat flux into the thermal shield) and the 
conduction response of the eroding solid. 

To illustrate the effects of the two variables con- 
sidered in this paper, boundary motion and variable 
conductivity, we arbitrarily choose a standard re-entry 
trajectory and a material which will erode when sub- 
jected to a high-temperature, high-velocity gas flow. 
Graphite is characteristic of materials which experi- 
ence mass loss by combustion and sublimation in such 


Application to Re-Entry Problems 
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Fic. 1. Re-entry body performance curves. 


environments. In addition, the physical and chemical 
data required for numerical calculation are readily 
available. Furthermore, the variation of graphite 
thermal conductivity with temperature is sufficient to 
produce a significant change in the temperature pro- 
file from that calculated using a constant value for 
conductivity. 

For the heat-flux history at the wall, the convective 
heat-transfer rate at the stagnation point of a blunt- 
nosed hypersonic re-entry body is used. This heat- 
flux schedule was originally outlined in reference 11 
and presented in terms of dimensionless time in refer- 
ence 12. 
material, the relations of reference 13 are adopted. 
The heat flux into the body is written: 


To compute the rate of removal of surface 


ee Ue(Cy. /2) | hse + Clo. Oree P= | — 


[tw ‘Pell, (Cp 2) )A po} 25) 


where awIo-Qpow is a measure of the enthalpy increase 


~ OUTER WALL 









__ 
4000 + My * 
\ 
3000 + 
DENISON & 

. DOOLEY 
ie — — — PRESENT 
4 ANALYSIS 
2000+ 
e K= 3360/T9-6 
a 
= 
ey 
S 

1000 + 

INNER WALL 








ie) 5 10 15 20 
TIME, SECONDS 


introduced due to simplification of heat flux 
computation. 


Fic. 2. Error 


SCIENCES—JANUARY, 1960 

due to combustion, and [My/p.t.(Cye/2) |Ayo iS a Meas- 
ure of the energy absorbed in the vaporization process. 
For comparative purposes, we may assume that these 
two terms are equal. The term /,,, is the enthalpy of 
air at wall conditions of temperature and dissociation 
and is small compared to the total enthalpy h,,. A 
close approximation for large Mach Numbers is 


(se — New) = (ur?/2gJ) 


and, with this simplification, the heat flux may be 
written: 


d = Pdlle (Cye/2) (ur? /2gJ) 


where p,u,(Cy/2) is a characteristic flow parameter. 
The terms p,, “,, and cy are the gas density, gas ve- 
locity, and local skin friction coefficient, respectively, 
based on conditions at the outer edge of the boundary 
layer. For trajectories of current interest in the ballis- 
tic missile field, m,./p,u.(Cz-/2) is approximately a con- 
stant* over a large temperature range and, for graphite, 
is equal to0.17. The mass flux at the wall is then given 
byt 
My = 4250 g/ (7/2) 
and the boundary motion is prescribed by 
(26) 


ds/dt = My/p = 40.4 g/(m7/2) 


The following parameters are chosen for calculating 
an illustrative re-entry trajectory. 


Up = 20,000 ft./sec., re-entry velocity 
CpA/m = 0.02 ft.?/slug, drag/ weight ratio 
On = 40°, re-entry angle 

R = | ft., nose radius 


The resulting heat flux, free-stream velocity, and 
boundary motion as a function of time from 200,060 ft. 
to impact are shown in Fig. 1. For all temperature 
profiles, a zero heat-flux boundary condition at the 
inner surface and a 0.0833-ft. wall thickness are as- 
sumed. 

To estimate the magnitude of error introduced by 
the restrictive assumptions used to calculate heat flux 
and erosion rate, comparison is made in Fig. 2 to the 
laminar aerothermochemical analysis of Denison and 
Dooley.'* The latter follows Lees’ stagnation-point 
theory but includes the effects of combustion and dis- 
sociation. Since the same heat-conduction solution 
is used for both curves shown, the difference betweer 
the two is due to the assumptions introduced betweea 
Eqs. (25) and (26) and the assumptions used in refer- 


* This assumption limits the analysis to surface combustion 
where the rate of mass removal is primarily controlled by the 
amount of oxygen at the surface and implies 0 and 0, accommo- 
dation coefficients of unity. For a more detailed statement of 
the nature of the combustion-sublimation process, including 
cases where the accommodation coefficients are less than unity 
and the rate of mass removal is determined by the sublimation 
rate, see reference 14. 

+ The authors are indebted to M. R. Denison of Aeronutronic 
for suggesting this relation. 
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TEMPERATURE PROFILES 


ence 11. Close agreement is obtained until the outer- 
wall temperature reaches 3000°R.; neglecting the 
effects of wall temperature and mass flux on the con- 
vective heat transfer rate beyond this temperature 
overestimates the heat flux and yields an outer-wall 
temperature approximately 10 per cent higher than 
the more exact solution shown. The inner-wall tem- 
peratures agree quite closely due to the short re-entry 
time. For cases with smaller re-entry angles, the agree- 
ment would be less favorable. 

The main comparisons to be made in this section in- 
volve the effects of conductivity variation and erosion 
rate. Fig. 3 shows the outer- and inner-wall temper- 
ature as computed for various conductivity relations. 
The wall thickness used was 0.0833 ft., and the vari- 
able conductivity for graphite was assumed to follow 
the relation K = 3360/(T)°*, a good approximation 
between 1500 and 6000°R. For the given heat and 
mass-flux schedules, a value K = 30.0 gives the same 
maximum outer-wall temperature as that computed for 
a variable K. The differences between the two his- 
tories may be accounted for by the fact that the con- 
stant conductivity initially underestimates thermal 
diffusion into the solid, but, during the final part of the 
trajectory, the constant conductivity is too high and 
the rate of temperature rise at the outer surface is de- 
creased. Corollary reasoning may be used to compare 
the two inner-wall temperature schedules. For equal 
maximum temperatures with this trajectory, the con- 
stant conductivity is only slightly less than that corre- 
sponding to the mean outer-wall temperature, indi- 
cating that the influence of the outer surface on con- 
duction is dominant. 

The effect of boundary motion on the temperature 
profiles is shown in Fig. 4. For cases such as the 
present where time is short and boundary motion is 
small (less than 10 per cent of the original thickness), 
any differences in the computed inner-wall temper- 
ature are well within the inaccuracies of the conduc- 
tivity data. For materials such as silica, where the 
mass loss represents a significant fraction of the original 
thickness, only a low conductivity would keep the rear 
wall temperature from being affected. 

Outer-wall temperatures decrease as the erosion rate 
increases. At any instant of time, the exposed surface 
is at a lower temperature than the surface which was 
removed. In cases where erosion is severe and con- 
ductivity high, the decrease in effective thermal ca- 
pacitance will tend to increase the outer-wall temper- 
ature over that obtained with no erosion. When mate- 
rials are used which show a strong increase in mass 
flux with temperature, this process tends to be self- 
limiting. It has been well documented (cf. reference 
16) that injection (or, equivalently, boundary motion) 
decreases the local Stanton Number. Thus, any large 
increase in surface temperature would decrease the 
heat flux to the wall and limit the surface temperature 
rise. 

Finally, mention should be made of the effects of wall 


curvature. In cases such as the present where the 
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Fic. 3. 


thickness /diameter ratio is small, effects on the tem- 
perature distribution are of lower order. The radial 
case is presented here as the more general solution, 
the flat plate solution requiring only a simple change in 
the coefficients 8;, y;, and yw; appearing in Eqs. (17)- 
(19). 
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Authors’ Note 


Subsequent to the submission of this paper for publi- 
cation, a solution to the one-dimensional melting slab 
problem with constant conductivity was presented by 
Citron.'® In this reference, the boundary condition at 
the moving face is constant temperature, with a zero 
heat-flux condition at the rear surface. Citron’s semi- 
analytical solution for arbitrary heat input at the ex- 
posed face solves explicitly for the boundary motion by 
requiring that the given heat flux Q(t) at the exposed 
face supply both the heat of fusion of the melted mate- 
rial, and the heat absorbed by the solid. Although 
analytical in its formulation, Citron’s analysis requires 
the numerical solution of a complicated and highly 
nonlinear ordinary first-order differential equation in- 
volving the location of the moving boundary and its 
time derivative. The numerical solution of this equa- 
tion is of the same order of complexity as the complete 
numerical solution to the transient conduction prob- 
lem as given in the present paper. To consider the 
case of finite latent heat (or, more specifically, the case 
in which the difference between the heat of ablation and 
the heat of reaction of the ablated material is not 
equal to zero) in the present solution, g, is found at 
each time point by subtracting the net required heat of 
ablation from the total heat flux to the wall. Insofar 
as the required heat of ablation depends upon boundary 
motion, it is again determined by averaging over the 
time interval using the previous step in the iterative 
solution of the difference equation. 
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Comments on ‘‘Limiting Circulatory Lift of a 
Wing of Finite Aspect Ratio”’ 


G. J. Hancock 
Dept. of Aeronautical Engineering, Queen Mary College, 
University of London, London, England 


July 6, 1959 

T° REFERENCE 1, McCormick indicates the existence of the 
maximum circulatory lift of a finite wing due to its trailing 

There are two points worth further discussion. 

accuracy of McCormick’s formula that 


vortex system 
The first concerns the 


(Cr)maz = 1.21 (aspect ratio) (1) 


The second point questions the direct application of Eq. (1) to the 
finite jet-flapped wing 

In the derivation of Eq. (1), 
elliptic spanwise loading distribution and that the trailing vor- 
ticity lies in a plane inclined to the free stream at the induced 
downwash But, the 
ultimate downwash angle far behind the wing will be twice the 
Thus, the effects of the curvature of the 


it is assumed that the wing has an 


angle calevlated at the wing lifting line 
inclination at the wing 
trailing vorticity in the stream direction are neglected 
it is possible to include these effects by considering the forces on a 


However, 


control surface at infinity in the usual manner.? A wing, which 
has an elliptic load distribution given by 
1/2 (2? 


r = Mol — y2/s?)'? 2 


together with its associated trailing vorticity, has the force dis- 


tribution shown in Fig. 1 





The lift coefficient, defined as the force derivative normal to 
the free stream, is given by 
Ci tA7(1 — 27? (3 
where A aspect ratio and 4 I)/4sV. The maximum lift 
coefficient occurs when y = 1/+/6 so that 
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(1) Five, double-spaced, typewritten, manu- 
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(Cx) maz = 0.855 A (4) 


This formula gives a value 30 per cent less than that of Eq. (1). 

The second point mentioned above is the application of Eq. 
(4) to the case of the maximum circulatory lift of a finite jet- 
flapped wing, which was established experimentally by Lowry 
and Vogler,’ since the downwash at infinity downstream of the 
wing depends on the total lift of the wing-jet system. The finite 
jet-flapped wing, assuming elliptic spanwise loading on both wing 
and jet, with its associated trailing vorticity is shown in Fig. 2. 

In Fig. 2, lw represents the wing-bound vorticity, I'y repre- 
sents the jet-bound vorticity due to the streamwise curvature, 
and J is the jet thrust. Hence, the lift coefficient is given by 


Cr = rAy7(1 — 2y?) + 2yCy (5) 


where y = (Tw + I'y)/4sV and Cy = J/(1/2)pV?2S. 
tory lift coefficient is now defined as 


Crp = Ci — Cysin ¢ (6) 


The circula- 


where ¢ is the angle of the jet to the free stream. Substituting 
Eq. (5) into Eq. (6), the maximum value of Cup can be calcu- 
lated by considering the variation of y, keeping Cy and ¢ con- 
stant. Hence 
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(Cry) = 0.855A [(1 + 2C,/rA)?!2 — 1.83(2Cy/xA) sin ¢] 


mar 


This equation only gives the maximum circulatory lift possible 
for a given value of C,; it does not state that this is attainable, 
However, in a range of experiments? of a series of wings (A = 
2.8, 5.6, and 8.4) and thrust coefficients (upper limit of Cy about 
7A/2) with sin g = 1, maximum values of Cip were found. It 
is interesting to note that the variation of C; between 0 < Cy; < 
aA /2 in this case does not materially change the overall values 
in Eq. (7). The comparison of experiment and the theory of 
Eq. (7) is shown in Fig. 3. 

The theory underestimates the experimental values in all 
cases, which may be due to two causes: 

(1) Because of suctions at the spanwise edges of the jet, the 
span of the jet together with the trailing vorticity, will increase, 
so increasing the effective aspect ratio. This effect will be im 
portant for the low aspect ratios. 

(2) If the loading changes from the elliptic distribution, the 
average downwash at infinity will decrease, so increasing ( Ciy Pa 


This effect will be important for the high aspect ratios. 
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Author’s Reply 


Barnes W. McCormick, Jr. 


Associate Professor of Aeronautical caplmeciog, 
The Pennsylvania State University, University Park, Pa. 


August 27, 1959 


T IS GRATIFYING to see that my note has stimulated additional 
work along these lines. Dr. Hancock has appeared to re- 
solve very neatly an assumption in my original work and has 
presented an interesting extension of the theory for the jet 
flapped wing. 

Another possible reason for the refined theory underestimating 
the experimental results can be obtained from the fact that the 
approximate solution agrees as well as it does with experiment 
(Compare Eq. (1) and Fig. 3 of Hancock’s paper). 
the effect of viscosity and the inherent instability of the trailing 
vortex sheet are such that that portion of the sheet adjacent to 
the wing would contribute relatively more to the induced velocity 
than would be obtained assuming a sheet of infinite extent. 


Certainly, 
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On the Maximum Attainable Temperature 
in Resonance Tubes 


Ascher H. Shapiro 

Professor of Mechanical Engineering, Massachusetts Institute o 
Technology, Cambridge, Massachusetts 

July 16, 1959 


~ 


SYMBOLS 
r = speed of sound 
f = frequency of oscillation 
k = ratio of specific heats 
L = length of resonance tube 
p = pressure 
T = temperature 
t = time 
u = particle velocity 


Vw = wave speed 





(7) 
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x = longitudinal distance 
p = density 


INTRODUCTION 


ere AND RESLER! suggest that the maximum temperature 
of the air in a resonance tube is limited by the condition 
that a continuous pressure wave entering the inlet barely steep- 
ens to a shock wave when it leaves at the exit. At this condi- 
tion, viscous dissipation within the tube ceases, and, apart 
from frictional heating in the boundary layer, there is no longer 
energy input to the tube. 

The analysis of reference 1 is based on the assumption that the 
temperature in the tube is uniform along its length. Experi- 
mental data? show, however, that there is a roughly linear 
temperature gradient within the tube. This suggests an addi- 
tional mechanism for inhibiting energy dissipation within the 
tube, namely, that a time-mean temperature gradient can pre- 
vent a wave from steepening. This mechanism is analyzed be- 
low and leads to an estimate of the maximum attainable tem- 


perature at the closed end. 


ANALYSIS 

Let there be a basic, time-invariant, longitudinal temperature 
distribution in the tube, 7(x), shown by the heavy line in Fig. 
l(a). For purposes of analysis, this is replaced by the dashed, 
step-like curve involving finite temperature jumps AT. Note 
that any unsteady temperature fluctuations due to the passage of 
waves are superposed on 7(x). 

When an incident wavelet 9 impinges on a temperature dis- 
continuity, there results a transmitted wave 3 and a reflected 
The resonance tube is, of course, filled with a 


wave &, Fig. 1(b). 
Moreover, there is 


complex pattern of waves of both families. 
no possibility of ‘simple waves” in the presence of the longitu- 
dinal gradient 7(x). Nonetheless, it is instructive, and quite 
likely indicative, to analyze a succession of waves of the type of 
Fig. 1(b) from the point of view of their steepening, without tak- 
ing account of the leftward-moving system of waves (except for 
the waves & reflected from the temperature discontinuities). 
The dynamical equations for the three wavelets are: 


p2 — pi = pici(te — 4) (1) 
ps — po = —prl(uz — U2) (2) 
Ps — ps = patsy — U4) (3) 


while the matching conditions across the temperature discontin- 
uity are: 
ps = Ps 


UW = Uy, pr = Ps, Us = Us, 


The 9 wave is assumed to be of small amplitude, and, thus, 
3. and ® are also small. Moreover, 7(x) is so divided that the 
amplitude A7/T is of the same order as the wave amplitude 
(ps — pi)/pi. Then, ignoring quadratic terms in these ampli- 
tudes (compared with unity), one may solve Eqs. (1)-(3) and 
arrive at 
(us; — U2) (Ue — uu) = —( pss — pici) 2p) (4) 
Now, for a perfect gas, the acoustic impedance is given by 


d(pc)/pc = (dp/p) — ((1/2)(dT/T)] 


Therefore, Eq. (4) becomes 


(u3 — U2)/(u2 — m) = (Ts — T1)/4T1 (5a) 
and it then follows easily from Eqs. (1), (2), and (3) that 
(ps — p2)/(p2 — fr) = —(Ts — T1)/4T) (5b) 
(us; — us)/(ue — wm) = 1 + (Ts — T1)/4T%) (5c) 
(ps — Ps)/(p2 — fr) = 1 — (Ts — T1)/4T) (5d) 


If Ap is the pressure jump across the forward-moving wave as 
it traverses the gradient 7(x), Eq. (5d) shows that it follows the 
law (Ap)'T = constant. 

Consider now a second wavelet propagating rightward. Does 
it overtake or fall behind the first wavelet? If 5V,, is the excess 
of propagation speed of the first wavelet over the second, reckoned 
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at a fixed time, then 
—5Vy = (uz + C3) — (Ug + 4) 


Now it may be shown that (u3 + c;) differs from (#2 + ce) only in 


second order. Making use of the previous results, we then find 


ane (6V, a) = \(k + 1) 2] [( tte — 4) C1] = 
(1/2) ((T% — 1T%)/T1) 
or, in the limit, 
5V,,/ix = |[(k + 1)/2](Ou/dx) + (c/2T)(dT/dx) (6) 
For 6V,,/6x to be zero, thus inhibiting further wave steepening 
and viscous dissipation, the required temperature gradient is, 
therefore, 
(1/T)(dT/dx) = —[(k + 1)/c](Ou/dx) (7) 


ORDER-OF-MAGNITUDE CALCULATIONS 
Proceeding now with some order-of-magnitude estimates, we 
may write, for a progressive wave, 
Ou/Oox = (1/c)\Ou/Ot 


Then, since the tube resonates at approximately quarter wave- 
length acoustic frequency, and since, for strong resonance, the 
particle velocities are of the order of the jet velocity and the lat- 
ter is of the order of the speed of sound, the maximum 0u/Ox is of 
order 

(0u/OX) mar ZS (1/c)(0u/O) nar LS 8f S 2c/L 


Thus, finally, 
AT/T 2 Ak + 1) (8) 


This means that in the absence of energy loss through the tube 
walls, the maximum temperature rise from open to closed end 
would be several times the mean temperature. The maximum 
rise so far obtained in experiments is about twice the mean 
temperature. 

These results are evidently only indicative of orders of magni- 
tude. Indeed, if one repeats the calculation for the condition 
that the gradient 0u/dx shall not increase while following the 
wave, it is found that the required d7/dx is just twice that given 
by Eq. (7). 

It seems likely that the actual limiting mechanism is a combina- 
tion of the two factors already mentioned: (1) as the average 
temperature in the tube increases, waves require a longer dis- 
tance to steepen, and (2) a longitudinal temperature gradient it 


self inhibits wave steepening. 


REFERENCES 
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461-2, July, 1959 
2 Sprenger, H., Uber Thermische Effeckte in Resonanzrohren, Mitteilung aus 
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Solution to the Flow About a Two-Dimensional 
Flat Plate at Infinite Mach Number? 


Donald L. Chubb* 
Princeton University, Princeton, N.J. 
July 17, 1959 


— OBTAIN A SOLUTION to the flow about a flat plate aligned 
perpendicular to the free stream at infinite Mach Number, 
a mathematical technique known as the ‘‘method of integral 
relations” was used. This method was first proposed by Dorod- 
nitsyn; a thorough explanation of Dorodnitsyn’s method is 
given by Hayes and Probstein.! In the solution, it has been 
assumed that the flow is homoenergetic and that the gas is 
calorically perfect. 

The solution obtained for the flat plate, using the method 
of integral relations, is for the subsonic region between the 
shock and the body only. A solution cannot be obtained for 
any point beyond the edge of tt:z plate, since at the edge there 
is a discontinuity in the flow. This discontinuity is a sonic 
point, and the velocity derivative becomes infinite as a result. 
If a subsonic region occurs outside the edge of the plate, it will 
not be taken into account; however, since the region is subsonic, 
it will influence upstream conditions. Therefore, as a result of 
not including this region, the solution will be incomplete. 

An infinite Mach Number was chosen for the solution in order 
to simplify the problem. As a result, the solution obtained 
will be true for the Mach Number regime where the Mach 
Number independence principle is applicable. In taking the 
limit 1/,, ~ o, the free-stream velocity, U’,, and density, p., 
remain finite; and the free-stream temperature, 7’, 
sure, p,,, approach zero. 


, and pres- 


+ This article resulted from independent work done during the author's 
senior year at Princeton, under the direction of Professor Wallace D. Hayes 
* The author wishes to acknowledge the helpful discussions with Richard 
solution 


Sebass and the assistance of Joseph Burns in programing the 


for the IBM 650. 
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In the method of integral relations, the region between the 
shock and the body is divided into N equally spaced strips. 
For the solution of this problem, N=1 was chosen. Therefore, 
there is only one strip with its boundaries being the shock and 
the plate. More exact results could be obtained if a larger 
value of N were chosen. Belotserkovskii,? using the method of 
integral relations for the flow about a circular cylinder, found 
that the results for VN = 2 and N = 3 agree closely. 

In this case, the partial differential equations that govern the 
nature of the problem may be reduced to total differential 
equations, if several functions are represented by interpolation 
polynomials of the form 

\ 
P,(x, ¥, 4 ...Un) = to an (x)y™ 
m=0 
where the a,,(x) depend linearly on the values of the functions 
P, on the strip boundaries and the u terms are dependent func- 
tions of x and y, such as velocity and pressure. It is obvious 
that, for our case, where N = 1, these polynomials will be linear. 
Thus, in the solution we have made the approximation that 
several functions vary linearly between the body and the shock. 

As a result of defining these interpolation polynomials, we 
total differential 
a numerical integration, the problem can be 


have three equations in three unknowns. 
Therefore, by 
solved. By using a finite-difference process, the solution was 
programmed for an IBM 650. The results of this solution for 
1.4 are plotted in Figs. 1-3. 


standoff distance was found to be Ay = 


The nondimensional 
0.8030. This 
dimeasional standoff distance is defined as Ap/Xo, where Ap is 


air with 7 = 
non- 


the distance between the shock and the body at X¥ =0 and Xo 
is half the width of the body. 

For this solution, where N = 1, the variables of the problem 
such as the velocity distribution can only be determined along 
the plate surface and immediately behind the shock. Thus, 
the sonic line is determined by only two points and, therefore, 
is a straight line as is shown in Fig. 1. Since results can be 
obtained oaly up to the edge of the plate, the location of the 
sonic point on the shock had to be determined graphically by 
continuing the curve in Fig. 2 beyond the edge until sonic velocity 
occurred. 

It can be seen from Fig. 1 that a small subsonic region occured 
outside the edge of the plate. This indicates the solution is 
incomplete; but since the velocities in the region are nearly 
sonic, and the region is small, as well, it should have only a 
slight effect upon the results. 


REFERENCES 
1 Hayes, Wallace D., and Probstein, Ronald F., Hypersonic Flow Theory, 
Academic Press, 1959. 
2 Belotserkovskii, O. M., Flow Past a Circular Cylinder With a Detached 
Shock, Vychisl. Mat., 1958, No. 3, pp. 149-185 
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Boundary-Layer Stabilization by Distributed 
Damping* 


Max O. Kramer 
Coleman-Kramer, Inc., Los Angeles, Calif. 
May 26, 1959 


I JuNE, 1957, test results were reported in the JOURNAL which 
indicated that a resilient energy-absorbing coating can delay 
the boundary-layer transition from laminar to turbulent. Since 
these preliminary tests had been conducted under unusual flow 
conditions, they were later repeated under more classical condi- 
tions, such as stepwise increased speed and zero pressure gradient 
in all directions. 

Fig. 1 shows the dimensions of the ducted rubber coating that 
It consisted of a heavy diaphragm sup- 
The space between 


was tested recently. 
ported by a multitude of tiny rubber stubs 
the rubber stubs could be filled with damping fluids of various 
viscosities. Various types of this coating were tested. These 
types differed in their elasticity module or stiffness as well as in 
their inherent damping. The stiffness of the various coatings was 
measured in Ibs. per cu. in., meaning the load on 1 sq. in. of the 
coated surface that would lead to a compression of 1 in. of a rigidly 
supported flat sample of the coating if the coating behaved 
linearly up to such large deflections. The actual deflections 
during such stiffness tests were in the range of 5/1,000 in. to 
10/1,000 in. The inherent damping of the coating was deter- 
mined by the percentage loss in amplitude after one bounce when 
dropping a cylindrical piston of 1/2-in. diameter and 1/2-o0z. 
weight from a 2-in. height onto the surface of a rigidly supported 
flat sample of the air-filled coating. The damping values men- 
tioned below are valid at 75° F. 

The coated body tested consisted of a cylindrical section 26.5 
in. in length and 2.5 in. in outside diameter, headed by an 18.5-in. 
tip of the tabulated contour. The first six inches of the tip, 0.5 in. 
at the transition from the tip to the cylinder and 1.5 in. at the 
aft end of the cylinder were not coated. The coated models were 
sprayed with airplane dope and polished in order to give them a 
high-gloss finish. 

The coated body or model was mounted on a sting in front of a 
cylindrical body of 2.5-in. diameter. The model plus its support- 
ing body were towed in water up to a speed of 30 knots. The 
ambient turbulence was smaller than 0.1 per cent. A 
strain-gage spring permitted the measurement of the drag of the 
model independently of that of the supporting body. All meas- 
urements were made as comparative measurements between the 
coated models and an uncoated model of equal shape having a 
high-gloss rigid surface. All results were evaluated in the con- 
ventional way as drag coefficient over Reynolds Number, where 
the drag coefficient was referred to the outside surface of the 
model and the Reynolds Number to the length of the model. 

Fig. 2 presents the results obtained. The ‘‘fully turbulent” 
and ‘fully laminar’ curves are the classical results on surface 
friction—i.e., the Shoenherr curve and the Blasius law 

Curve A refers to the rigid reference model with its high-gloss 
rigid surface. The performance of this model was normal, being 
almost fully turbulent at the maximum Reynolds Number, R = 
15 X 108. 

Curve B refers to a coated model with 1,600 Ibs 
stiffness and 44 per cent inherent damping of its coating. In 
reference to the rigid model, Model B reached 35 per cent drag re- 
duction at R = 15 X 108. The viscosity of the damping fluid for 
best performance was 1,200 centistoke. 

Curve C refers to a coated model with 800 Ibs. per cu.in. stiff- 


per cu.in. 





* Complete details of this work will appear in the February, 1960, issue 


of the Journal of the American Society of Naval Engineers 
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(dimensions in 1/1,000 of an in.). 


ness and 47 per cent inherent damping of its coating. In ref- 
erence to the rigid model, Model B reached 59 per cent drag re- 
duction at R = 15 X 108 The viscosity of the damping fluid for 
best performance was 300 centistoke 

Curve D refers to a coated model with 600 Ibs. per cu.in. stiff- 
ness and 57 per cent inherent damping of its coating. In reference 
to the rigid model, Model D reached 10 per cent drag reduction 
at R = 15 X 108 Within the range from 10 to 200 centisoke, 
the viscosity of the damping fluid had no effect on the per- 


formance of this model. 
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Fic. 2. The drag coefficient of various models as a function of 
the Reynolds Number. Curve A is the rigid reference model, and 
Curves B, C, and D are fully coated models with a stiffness of the 
coating on the cylindrical section ( B = 1,600 Ibs. per cu.in., C = 
800 Ibs. per cu.in., and D = 600 Ibs. per cu.in 
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Analysis of Further Data on the Effect 
of Isolated Roughness on Boundary-Layer 
Transition in Supersonic Flow 


Darwin W. Clutter* and A. M. O. Smith** 
Douglas Aircraft Company, Inc., El Segundo, Calif. 
July 20, 1959 


is REFERENCE I, the authors proposed a new physical model to 
describe the effect of isolated roughness on boundary-layer 
transition in compressible flow. At the time reference 1 was 
prepared, only a limited amount of data was available for check- 
ing the correlation presented there. A recent publication by 
Van Driest and McCauley? supplies further data on the effect 
of three-dimensional roughness elements on transition, and 
these data further support the correlation of reference 1. 


It was shown in reference 1 that a roughness element would 
cause premature transition if its local Reynolds Number exceeded 
a certain value. This value is called the critical roughness 
Reynolds Number. In subsonic flow, the roughness Reynolds 


Number is defined as 


Rx = upk/ vy, (1 ) 
wherek = height of roughness element 
u; = velocity at a height k in the boundary layer in the 
absence of the roughness 
vy, = kinematic viscosity at height k in the boundary layer 


The value of the critical roughness Reynolds Number was shown 

to be a function of Mach Number and roughness shape only. 
When the local flow about the roughness is supersonic, the local 

Reynolds Number presented in reference 1 is defined as 


Ry* = [0(T*, ps)a*k] /[u(T*)] (2) 


where p(7™*, ps) is the density based on 7*, the sonic temperature 
of the stream filament at the roughness height, and on p35, the 
local pressure outside the boundary layer; a* is the sonic ve- 
locity at the height &, and yu is the viscosity at temperature 7*. 
The argument for sonic conditions is that for a blunt body in 
supersonic flow, sonic conditions occur near the top of the rough- 
ness. The density is based on the pressure outside the boundary 
layer because experiments show no great change in static pres- 
sure near a roughness element. The critical values of R,.* seem 
to be nearly constant for a given shape of roughness at all super- 
sonic Mach Numbers. 

Values of critical roughness Reynolds Number are presented 
in Fig. 1. The values of R,.1/2 or R,.*!/? are plotted against the 
Mach Number at the top of the roughness J, = u,;/a;. The 


* Aerodynamicist, Aerodynamic Research, El Segundo Division. 
** Supervisor, Aerodynamic Research, El Segundo Division. 





Variation of critical roughness Reynolds Number Rx 
and Ax* with local Mach Number. 
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Fic. 2. The roughness Reynolds Number required to produce 


transition at the roughness, as a function of the local rough- 
ness Mach Number. 


recent data of Van Driest and McCauley (for spheres) give 
nearly the same values of critical roughness Reynolds Number 
as were found earlier for two-dimensional roughness (wires) by 
Van Driest and Boison.* The diagonally shaded rectangle and 
the dashed line are subsonic data from reference 1. As defined 
above, the velocity u, in the formula for R, is the velocity at a 
height & in the boundary layer in the absence of the roughness. 
One would expect the true velocity at the top of the roughness 
element to be about twice this value. The rectangle of dashed 
lines represents values of low-speed R,’s based on the ‘‘true ve- 
locity.’’ As is to be expected, these R,’s are more in line with 
R;.*, which is also based on the velocity near the top of the rough- 


ness. 


Fig. 2 presents similar values of the roughness Reynolds Number 
required to produce transition at the roughness. The trend of 
the data is the same as that of the data for critical roughness 
Reynolds Number, but, as is to be expected, the Reynolds Num- 
ber required to produce transition at the roughness is larger 
than the critical value. The single point in the diagonally shaded 
region is the roughness Reynolds Number found for spheres in 
low-speed flow by Klebanoff, et al. 


It should be noted that most of the supersonic data presented 
here on roughness have been obtained in the Jet Propulsion 
Laboratory 12-in. or 20-in. wind tunnel. Furthermore, it should 
be noted that for a given Mach Number the transition Reynolds 
Number observed in the JPL tunnels remains essentially constant 
for varying unit Reynolds Number. Recent publications 
of transition Reynolds Numbers in other supersonic tunnels show 
rather strong variations with unit Reynolds Number (for ex- 
ample, reference 4). This variation may be due to tunnel 
turbulence. The data on transition due to roughness presented 
in reference 4 were also examined but were not sufficient to allow 
the roughness Reynolds Numbers as defined by Eqs. (1) and (2) 
to be determined either for critical roughness or for the case 
when transition occurs at the roughness. It would be of interest 
to see if roughness data from supersonic tunnels other than those 
of JPL also fit the correlation presented in reference 1. 
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An Experiment on Oscillatory Shock-Wave 
Boundary-Layer Interaction 


Wladyslaw Fiszdon* and Erik Mollo-Christensen** 
Massachusetts Institute of Technology, Cambridge, Mass. 
July 22, 1959 


INTRODUCTION 


pointed out the possibility of self-excited 


oe has 
oscillations of a boundary-layer interaction. 


The theoretical considerations which led to this result involve 


shock-wave 


rather radical approximations and simplifications which need 
experimental verification. 

We have recently performed some preliminary 
which appear to confirm the existence of a resonant frequency 
of a shock-detachment region.?, The experiments neither confirm 
nor invalidate the details of the theory. 


experiments 


DESCRIPTION OF EXPERIMENTS 

The experiments were performed in the M.I.T. 8 in. X 8 in, 
Mach Number 2 supersonic wind tunnel; Fig. 1 shows the arrange- 
ment. A flat plate, spanning the tunnel, was installed. It con- 
tained numerous pressure holes on one side and was also fitted 
out with a piezo-electric pressure transducer (B) and a traversing 
boundary-layer pitot tube. The shock generator consisted 
of two plates spanning the tunnel, with a movable plate in be- 
tween. A rectangular cavity of variable depth was thus formed. 
The bottom of the cavity (the front end of the moving plate) was 
provided with a piezo-electrical pressure transducer (A). 

The flow behind the detached shock wave formed upstream 
of the shock generator will under certain circumstances be un- 
steady, as could be observed from the output of transducer A. 

Fig. 2 gives the frequency of the signal from transducer A vs. 
cavity depth /. The change in frequency with / is probably due 
to changes in mean temperature in the cavity, the shorter 
cavities running hotter than the long ones, in which the dissipa- 
tion by skin friction is more important than the irreversible 


processes in the shock waves going back and forth. 
RESULTS OF THE EXPERIMENT 


The maximum ratios between the power spectral densities 
of the inputs at transducers A and B were computed from the 


* Professor, Warsaw Technical University, Guest of M.I.T. 


** Associate Professor, M.I.T. 
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measured spectra for several Reynolds Numbers and are shown 
in Fig. 3. 

It should be remarked that the signals were almost simply 
harmonic, and therefore the maxima were very sharp. An 
apparent tunnel resonance at about 11,000 eps was ignored in the 
analysis. 

The result of Trilling’s theory is shown in Fig. 3 


CONCLUSION 
The experiment indicates that there is a resonant frequency of 


shock-wave boundary-layer interaction 
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Approximate Solutions for Supersonic Flow 
Over Wedges and Cones} 


A. G. Hammitt* and K. R. A. Murthy** 
Princeton University, Princeton, N.J. 
July 22, 1959 


HE SOLUTIONS for oblique shocks and conical flows are well 
known. Unfortunately, it is usually necessary to tabulate 
the solutions since simple analytic relations are not available 
It is the purpose of this note to point out some simple approxi 
mate relations which can be used at large ./ for these two cases 


OBLIQUE SHOCKS 


It is not generally pointed out that for oblique shocks the rela- 
tion between the shock angle, @,, and turning angle, 0,, can be 
written exactly as a cubic in terms of wave angle 


A tan* 6, + Btan?@, + Ctané@, + D = 0 (1) 
where 
A=1+ [(7 — 1)/2] M2, C= 14 [(yv + 1)/2]M;? 
B = —(M,? — 1)-cot @., D = cot é 
This relation can be solved for 6, directly 
tan 6, = (1/3A){2(B? — 3AC)'/?- cos [¢ + (nx/3)] — B} 
and 
he ; 9ABC — 27A2D — 2B3 
oe 2(B? — 3AC)3/2 
n =0 gives the strong-shock solution, and n = 4 gives the weak- 
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Fic. 1. Variation of shock angle with Mach Number and wedge 
4 


angle for y = 


shock solution. The detachment angle, 6,4, is given by 


cot? Ona = (P 2) + V/(P2 


where 


) + [AC2M,? — 1)5]_— (3) 


[4(.Mi? — 1)3]} [18AC( AG? — 1) + 


27A? — (M,? — 1)?-C?] 
This cubic equation can be reduced to a quadratic for large M/ 


in two ways. For 6, small, it gives the well-known result 


6: = [(y + 1)/4]O0 + Wi Cy + 1)/4]6.}2 


if tan? 9, term is neglected. 
For 6, large, 


+ (1/M,2) (4) 


it gives 


(M,\? — 1) cot 6% 
tan 6, 


2{1 + [(y — 1)/2] 4} 
(M,2 — 1)? cot? 6, | + ey + MP 
= « o) 
HL + [Cy — 1/242 1 + [Gy 2).M? 


found by dividing the cubic by tan 6, and neglecting the last term. 
A more exact relation can be found by an iteration method using 
the solution from Eq. (5) as an approximate value for the D/tan 


6, term and solving the resulting quadratic. These various 


relations are shown in Fig. 1. For M/ > 5, good results for all 


angles can be obtained by using the appropriate approximate 


formula. 


ConricaL FLow 


The conical-flow problem can be treated by a method similar 
to that of reference 3, but not limited to small angles. The 
velocities between the shock and the body can be expanded by 
the Taylor series to give 


u/uUe = 1 — h? + (1/3) cot 6-:h? — ag-ht +... | (6 
/ ») 
v/ue = —2h + cot 6.h? — 4ash? +... ) 
where h = (0 — @,). 6 is the conical ray angle, and c refers to 


conditions on the cone. 
The boundary condition at the conical shock can be written 


—v;/Us Sin 0s:COS 0, = 


(vy — 1)/(y + 1)] sin? 6, + {2/1(y + 1)-n2]} (7) 


where s refers to conditions behind the shock. 
Combining these equations and retaining terms of order h?, 
gives 
[2 — (y +s) sin? 6,] -h?, + 
J} =0 (8) 


2 sin 20,.-h, — [(y — 1) 


sin? @. + (2/M,? 
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h, must be small so that the higher-order terms are negligible, 


but @. has no restrictions. The solution to the quadratic is 


— sin 20, 
h, = - oe 
[2 — (y + 5) sin? 9, 
sin 26, V* f(y — 1) sin 26. + (2/3%2)) 
y ‘a + 5) sin? 6-] § j | [2 — (y + 5) sin? 6,] f 
(9 


The correct sign is the one which gives the smaller positive value 
of h For sin? @. = 2/(7 + 5), the coefficient of the h?, term 
goes to zero, and the solution is 

(y — 1) sin? @ + (2/M,?) 


h, = : 0 
2 sin? @, (10) 


The detachment angle can be easily found from the quadratic 


(vy +1) — [C7 + 5)/M?] 


sin? 0.4 = : = 
4+ (7 + 5)(y7 — 1) 


(11) 


The shock angle and detachment angles as calculated by Eqs. (9 
~(11) are shown compared with the exact values in Fig. 2 
all other flow quantities can 


and surface Mach Num- 


Based on the same assumptions, 
be determined. The surface pressure 


ber, for instance, are given by the relations 


Ps 2) ? y-1 
— - M2 sin’? @. — , 4 
A ty et 1+¥1 


] » (12 
) 1+ ly - 2] 17,2 sin? @, { 
and 
, V2 cos? 6,(1 + 2h?, ; 
M,? = ecrreee es Z (13) 
1 + [(4 — 1)/2]M,)2(sin? 0, — 2h?,-cos? 0, 
to the same order of accuracy as the shock-wave angle 
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for comparable Foelsch-type nozzles. DISTANCE DOWNSTREAM OF THROAT (FT 
Fic. 3. Centerline Mach Number distributions. 
DISCUSSION 
The CWT operates supersonically (1.8 < M, < 1.8) using a 
symmetrical, flexible nozzle 8.5 ft. high and 20.5 ft. long, which 




















j produces about 7 ft. of testing length. The nozzle family is of the 
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point.*~* This allows the plate to match the theoretical aero- 
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produce a curvature discontinuity. The modification is defined = 
rd 6,p/6i, which describes a region of partial-strength waves in 3 CASE Ip OS 
the characteristic solution. When 6,p/6; reaches 1.0, the partial- oe oOo WA — THEORY FOR NEW FAMILY 
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flow so that 6;/»; 
servative in normal practice. 


0.7, which is not thought to be very con- 
However, with very minor experi- 
mental corrections, the flow had errors in Mach Number of less 
than +0.005. 

After these nozzles were put into operation, it was found that 
more testing length was needed. The leagth which can be gained 
depends directly on how large @,p/0; and 0;/y; can be made. 
The limit of 6;p/0; is 1.0, and the limit of 0;/»; can approach 1.0. 
To investigate the practical limits, the six A/, = 1.5 nozzles 


shown in Fig. 2 were computed and tested. 


Three cases, all 
starting at the present nozzle, were considered: Case I, y,;’’’ 
1.0; Case II, 0:p/0; = 0.5 with 6;/v; < 
1.0; and Case III, 0;/v; = 0.722 (as for the present nozzle) with 
6, p/0; 1.0. 
Number distributions for the extremes of the three cases with the 
first M, = Using the jack influence curves,®> the 
best corrected distributions were computed; they show that 
Cases I and II can be corrected to AJ = +0.003, while Case ITI 
cannot be significantly improved near the influence of the in- 


continuous with 0;/r; 
Fig. 3 compares experimental, uicorrected Mach- 


1.5 nozzle. 


flection point. 

Based on these tests at 1/7, = 1.5, the tentative conclusions for 
the Mach-Number range 1.3 to 1.8 are as follows: 

(1) The Mach-Number distributions can be corrected to AM = 
+0.003 when @,:p/0; = 0.7, and to AM = +0.01 when 6;p/6; = 
1.0. 

(2) The expansion parameter, 9;/v;, can approach 1.0 without 
appreciable overexpansion in the test region. 

Following this study, the CWT nozzle family was revised and 


calibrated. The new family is defined by the conditions that 


yi'"’ is continuous and @;p/0; = 0.7. As a consequence of this 
M, = 1.4 to 0.86 


The resulting nozzles are similar to Case I and give 


choice, 0;/v; for the family varies from 0.95 at 
M, = 1.8. 
the gains in testing length shown in Fig. 4. 
Numbers are within +0.008. 


The errors in Mach 
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Incompressible Wedge Flows of an Electrically 
Conducting Viscous Fluid in the Presence of a 
Magnetic Field; 


K. T. Yen* 
Rensselaer Polytechnic Institute, Troy, N.Y. 
July 23, 1959 


ew PURPOSE of this note is to discuss the two-dimensional 
flow of an electrically conducting viscous fluid past a wedge 
in the presence of a magnetic field. The governing differential 
equations and boundary conditions are given and analyzed. 

The differential equations for this problem are! 


q-Vq = —(1/po)VP + (ue/po(V xH)xH + vV2q (1) 
nV?H = q:‘VH — H-Vq (2) 


t+ This research was supported by the United States Air Force through 
the AFOSR of the ARDC, under Contract No. AF49(638)-23. 
tion in whole or in part is permitted for any purpose of the United States 
Covernment. 
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V-q=0, V-H=0 (3) 
where yn. = (1/ope) It is noted that, when E is zero, the current 


Jis 
J =VxH = on.(q x H) (4 
(The symbols in the above equations are as defined in reference 1 
If the following stream function y is introduced 
yy = ax’f(n), 1 = cyx4 (5 
then 


, 


u = Oy/oy = acx’ 4 (6) 
—v = Oy/Ox = ax’“(pf — anf’) (7) 


In the above relations, a and ¢ are prescribed constants. The 
indices p and gq are to be determined to yield similar solutions 
Similarly, a magnetic stream function may be introduced 


YW = Ax? F(n) (8 
H, = dW/dy = Acx?-9F’ (9) 
—Hy = oW/dox = Ax?~'(pF — qnF’) (10 


In order that the viscous terms be of the same order as the 


inertia terms, it is necessary, in accordance with Eq. (1), that 


p+qz=1. Then 
f'" + (a/cev) | — g)ff’ + (1 — 2a) — f’)] - 
(A*yu./acpov)((1 — g)FF”’ + (1 — 2q)( Fo? — F’)| =0 (11 
P = (1/2)} A2c?u.( Fo? — F’2) — a%c?}x?-4a (12 
where F,? is a constant. Eq. (12) indicates the presence of 


a transverse pressure gradient 


From the first equation of Eq. (2), it is found that, to the order 
of x1 74%, 
(nc/a)F’"’ = (1 — a\ fF" — fF) (13 
and from the second equation of Eq. (2), to the order of x ~*4, 
(nef a)F’ =(1— qQy(fF’ — f'F) (14) 
Thus, Eq. (13) is implied by Eq. (14) 
It is convenient to choose 
a= cpr, (Ay, acpyv) = m? 
so that 
u = Ugxraf’ 
H, = Hox? ~¢F’ 
m?* = (Fu, Uo? po) 


Eqs. (11) and (14) become 


f'’ +10 — Off" + (1 — 29q)1 — f’?)] - 
m?((1 — qg)FF’’ + (1 — 2q)(Fo? — F’?)] = 0 (15) 
yF" = (1 — g)\(fF’ — f'F) (16) 
where y = (7-/v). These two equations form a determinate 


system for the functions f and F. When the solution for F is 
found, P can be obtained from Eq. (12). 

From Eq. (4), it is found that the current is directed in a direc- 
tion perpendicular to the x, y plane and has the magnitude, to the 
order of x!~*4, 


J = —(A4ue/pov?)x!—34 FF’ (17) 

On the other hand, Eq. (4) also yields 
j = op(uHy — vH,) (18) 
Since u = v = Oat yn = 0, it follows that, for finite op,, (within 


the approximation as implied in Eq. (4)) the boundary condition 


F"' = 0, 7» = 0 (19) 


must hold. This is, of course, consistent with the velocity 
boundary conditions 
f=f' =0, »=0 (20) 
The additional velocity boundary condition is 
f'aol, y—> @ (21) 
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To solve the differential Eqs. (15) and (16), two boundary 
conditions for F are required. If it is required that the magnetic 


field is ‘‘transverse,’’ the boundary conditions may be taken as 
F=1, »=0 (22) 


P30. - gas (23) 


The constant Fo? in Eq. (15) should be taken as zero. For the 
case g = 0, the boundary conditions used in reference 1 are quite 
different, although they appear to represent approximately the 
same physical problem. (If the magnetic field is strictly trans- 
verse, the governing differential equations would be different.) 
Another set of physically reasonable boundary conditions for F is, 


for example, 


F = 0, 7 = 0 (24) 
F’—-FKk=1, »—> (25) 
Here, the magnetic field may be called “‘longitudinal.”’ It can 


be seen that the differential Eqs. (15) and (16), under the bound 


ary conditions (24) and (25), would admit the solution f = F, 
provided that y = 0 and the condition F’ = 0 at » = 0 is ad- 
mitted. No solution appears possible when m? = 1 

When g = 1/2, the problem reduces to the flow over a flat 


plate. For a constant transverse magnetic field, this problem 
has been studied by Rossow.? It is planned to solve the flat 
plate problem as formulated in this paper. Further work is under 


way and will be reported later. 
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Skin-Friction Calculation for Turbulent 
Boundary Layers in Adverse Pressure 
Distributions 


Earl M. Uram* 

Scientist, Research Department, Electric Boat Division 
General Dynamics Corporation, Groton, Conn. 

July 29, 1959 


A’ PRESENT, there are no satisfactory methods of calculating 
skin-friction coefficients for turbulent boundary layers in 
adverse pressure distributions without a priori knowledge of 
either the velocity distribution or both of its integral parameters, 
momeatum thickness (@) and displacement thickness (6*). 
When these quantities are known, the two available methods 
for such boundary layers are the Ludwieg-Tillmann formula 
(reference 1) and the law of the wall. These methods are dis- 
cussed in detail in reference 2, where it is concluded that calcu- 
lations based on the law of the wall yield more reliable results. 
It is also shown in reference 2 that methods exist which allow 
accurate determination of 6, but no consistently reliable method 
is available for the determination of 6*. Therefore, a skin- 
friction calculation method involving known or readily deter- 
mined quantities is needed. In the following development of 
such a method, skin-friction coefficients obtained from law of 
the wall calculations applied to experimental velocity distribu- 
tions are used, since no direct force measurements are available 
for boundary layers in adverse pressure distributions. 
It has been shown (reference 1) that the skin-friction coeffi- 
cient for boundary layers influenced by arbitrary pressure dis- 
tributions can be expressed as 


* Formerly Research Engineer, United Aircraft Corp., East Hartford, 
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MOMENTUM THICKNESS VELOCITY RATIO, 7 


Fic. 1. Correlation of cy with y for constant Re 


cr = cs(y™Re" (1 


where y is the velocity ratio, U’/U,, corresponding to the semi 
logarithmic portion of the velocity distribution or its extension at 
a distance from the surface equal to the momentum thickness 
On the basis of available experimental data, it is shown in Figs. 
1 and 2 that m = 2.04 and n = —0.1155, respectively. The 
data indicating a behavior different from the equation in Fig. 1 
are directly attributable to only two particular data sources 
Taking these values for m and n, the data can then be plotted 
as shown in Fig. 3 to yield the relationship 

cy = 15.0 K 1073 y?-04Rg ~0- 1155) (2) 
It was found that, for almost all flows, 7 can be correlated with 
the quantity 6/c; (with 6 expressed in inches) as shown in Fig. 4 
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7 1 by considering y as a function of H as shown in Fig. 5. Sub- 
y2 04 Rp ollss ? 

ss oe = se : mae = stituting the equation representing the correlation of Fig. 5 in 
Fic. 3. Skin-friction coefficient as a function of y?-04Rg~ 1155 8 1 I 8 8 


For greater accuracy, it is desirable to represent the two-dimen- 
sional data in Fig. 4 in two regions. The equations representing 
the correlations in Fig. 4 when substituted into Eq. (2) yield for 


conical flows, 
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Fic. 5. Correlation of y with H. 


Eq. (2) yields 
15.0 K 10 3flog(8.05 }{!.818)] 2.04 


c= (6 


Ro: 1155 


In most cases, Eq. (6) should yield better than +5 per cent 


accuracy in the value of c,. 


REFERENCES 
en uber die Wandshub- 
NACA TM 1285, 


1 Ludwieg, H., and Tillmann, W., Untersuchung 
spannung in Turbulenten Riebungsschicten 
1950 

2 Uram, E. M., Investigations of Incompressible Turbulent Boundary Layers 


Report SR-0949-1 


translated as 


United Aircraft Corporation Research Dept September, 


1959 
—+>——_ 


Erratum—‘‘A Theoretical Investigation 
of the Interaction Between Shock Waves 
and Boundary Layers’’* 


M. Honda 
Assistant Professor, Institute of High Speed Mechanics, 
Tohoku University, Sendai, Japan 


October 20, 1959 


ON PAGE 674, in Eq. (2.19), ao/a, should be a;/ao. 
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On the Normal Shock Wave Attached to the 
Curved Surface 


G. E. Gadd 

Aerodynamics Division, National Physical Laboratory, 
Teddington, Middlesex, England 

July 30, 1959 


A NOTE with the above title by Takeo Sakurai was published 
in the Readers’ Forum in July, 1959. It pointed out an 
error in Zierep’s work,! showing that his analysis in fact implies 
that, in the absence of singularities, normal shocks cannot occur 
on a convex curved surface in inviscid flow for upstream Mach 
Numbers in the range 1 to 1.662. However, it was suggested 
several years ago by Emmons? that at the surface immediately 
downstream of the shock there might be a singularity similar to 
that which occurs in low-speed fiow past a surface with a dis- 
continuity in curvature. Such a singularity has been investi- 
gated by Koppenfels.* The present author has investigated‘ 
the application of a singularity of this type to the normal shock 
case and has shown that a modification of Koppenfels’ solution 
appears to be applicable immediately downstream of a normal 








> 


shor 
pos 


July 


un 





per cent, 


y of Ec 


ith addi- 
tributed 


1 


eriment, 


to vield 


ore 


uni- 


eference 
». Sub- 


‘ig. 5 ir 


er cent 


andshub- 
M 1285, 


y Layers, 
»tember, 


1e 


ished 
it an 
plies 
yecur 
fach 
sted 
itely 
ir to 


dis- 


1 


. 


esti- $ 


ted! 
1¢ ck 
tion 
‘mal 








EADERS’ 


shock on a convex surface, making the occurrence of such a shock 


possible in the Mach Number range 1 to 1.662 
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A Note on the Explosion Solution of Sedov With 
Application to the Newtonian Theory of 
Unsteady Hypersonic Flow{ 


N. C. Freeman* 
Research Associate, Princeton University, Princeton, N.J. 
July 30, 1959 
(1) THE SEDOv RESULT 

A” EXACT analytical solution of the equations of inviscid com- 

pressible unsteady flow has been given by Sedov (reference 
to the solution may be made through Hayes and Probstein') 
This solution is the similarity solution for a constant-energy point 
In view of the recent work on problems of hypersonic 
near 1 


explosion. 
flow in the limiting form of the ratio of specific heats y 
(references 2 and 3), it would seem useful to analyze the Sedov 
solution in this limit and inquire what form such a solution would 
take. 
various + 
It may be shown that the convergence to the limit is nonuniform 
It is also not difficult to show that the non- 


Einbinder,® in a recent note, has examined the solution for 


but does not mention the interesting case of y = 1 
over the flow field. 
uniform behavior exhibited here is that which one would expect 
from the Newtonian formulation as derived in reference 3 

Using the notation of reference 2, the Sedov solution can be 
shown to be of the form 


t/t, = (to/te)* 1) 
v/t = (7,/r,)* 2 
p/ps = (1/2)(ro/rs)21 + [r./rs]?) (3) 
p/ps = (1/2)(1 + [ro/rs]?) (4) 
T/T, = (re/ro)? 5) 


in the limit as y — 1, where we have chosen for convenience the 


axisymmetric case (7 = 1) with uniform density in front of the 


shock wave (w = 0). The parameter V has been eliminated from 


the equations as given in reference 1. Here, we have written 
+1). It is clear that the limit y = 
cannot be taken in Eqs. (1) and (2) without further consideration 


The symbols p, v, 7, p, and r denote pressure, 


e= (7 — 1L)/ly lore = 0 
of the value of 75. 
fluid velocity, temperature, density, and distance from the point 
denotes values at the 
shock wave, and o denotes the Lagrangian variable. We see that 
from Eq. (1), for r, for order one, The 
whole flow is therefore concentrated at the shock wave r = ? 

The only particles which span the region between the shock unt 
the center are those which were very near the center originally— 
except in a 


of explosion, respectively. The suffix s 


we have r/r, > lase— 0. 


i.e., with small 7. In a similar manner, we see that, 


thin region near the shock, the density is small, the temperature 


very large, and the pressure has exactly half the value it has at the 


+ This research was sponsored by the Office of Scientific Research, USAF, 
under Contract AF 49(648)-465. 
* On leave from the National Physical Laboratory, Teddington, Middle- 


sex, England. 


FORUM 77 


shock wave. The situation in the limit, therefore, is that of_a 
thin, dense layer of fluid at the shock 
tenuous high-temperature region which extends to the center. 


The essential nonuniformity 


wave preceding a very 


This result was noted in reference 4. 
of the result in the limit is exhibited by Eq. (1). 


(2) THE NEWTONIAN-THEORY RESULT 


In developing the Newtonian theory for unsteady flow, it is 


shown in reference 3 that, for the problem of a cylindrical piston 


moving rapidly into a gas at rest, the narrow region between the 


and the piston could be treated as a moving 


strong shock wave 


piston of varying mass. This corresponds to the limiting case of 


y = 1. A simple application of Newton’s second law then gives 
p = (1/2r,)(d/dt)\(dr,/dt)(re2 — ro®)} = 
(dr,/dt)? + (1/2r,)(d?r,/dt?)(r,? ro?) (6) 
where the notation conforms with that used in the previous sec- 
tion. Here, r = ¢,(t) is the equation of the shock wave, and 
r, is the Lagrangian variable. Attempting to obtain a similar 
solution in the Sedov manner, we put 
r = t} - 
and obtain 
p = (1/8r.2)}1 + (” r.)?} (7) 


(4r,2)~! from the shock-wave relation, we 


Using the Newtonian result that in the 


By noting that p, = 
obtain Eq. (4) directly 
limit as y — 1 isothermal conditions prevail along particle paths, 
p/p = Ppo/P (8) 
where the suffix 0 denotes values at the time the particle crossed 
We deduce 
b/ Ps 
which is Eq. (3). 
To calculate the physical coordinate r, 
to use the continuity relation 


the shock wave. 


2)[ro/rlX1 + [ro/ro]?) (9) 


however, it is necessary 


prdr = pordro (10) 


or 


p)rodro (10a) 


(1/2)? = Sf (p- 


at constant 7;. po is the density in front of the shock wave 
though, that if we wish to extend 
0, this integral will diverge logarith- 


is not hard to find. 


It is immediately obvious, 
the limits of Eq. (10a) tor, = 
mically. The reason that it does diverge 
This occurs because we have used Eq. (8), 
rectly we should have used the equation for constant entropy along 


whereas more cor- 
particle paths 


Po? (11) 


p/p” = pi 


Using Eq. 11, we therefore obtain 


reef” Ee er 


It is now possible to proceed to the limit in this result, provided 


(1 —e)/(1 +e) 


ror, (12 


one retains the limiting form in ro 
Thus, 
(13) 


r2/r.2 = De (r./r.)  dir.j7,) 


And, 


hence, 


which is Eq. 1. 

We have thus shown that the assumption of Newtonian theory 
that breaks down in the derivation of the constant-energy simi- 
larity solution of Sedov is that of isothermal conditions along 
particle paths. The Newtonian theory is able, by modification of 
this assumption alone, to give the correct nonuniform dependence 


in the limit as y > 1. 
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Constant-Temperature Magneto-Gasdynamic 
Channel Flow 


Jack L. Kerrebrock* and Frank E. Marble** 
Daniel and Florence Guggenheim Jet Propulsion Center, California 
Institute of Technology, Pasadena, Calif. 


August 5, 1959 


INTRODUCTION 


- THE COURSE of investigating boundary-layer flow in con- 
tinuous plasma accelerators with crossed electric and mag- 
netic fields, it was found advantageous to have at hand simple 
closed-form solutions for the magneto-gasdynamic flow in the 
duct which could serve as free-stream conditions for the boundary 
layers. Nontrivial solutions of this sort are not available at 
present, and in fact, as in the work of Resler and Sears,! the 
variation of conditions along the flow axis must be obtained 
through numerical integration. 

Consequently, some simple solutions of magneto-gasdynamic 
channel flow were sought, possessing sufficient algebraic simplicity 
to serve as free-stream boundary conditions for analytic investiga- 
tions of the boundary layer in a physically reasonable accelerator. 
In particular, since the cooling of the accelerator tube is likely to 
be an important physical problem because of the high gas tem- 
peratures required to provide sufficient gaseous conductivity, 
channel flow with constant temperature appears interesting. 
Some simple algebraic solutions for the case of a constant tem- 
perature plasma are developed in the following paragraphs 


ANALYSIS 


If heat conduction and viscous effects are neglected, magneto- 
gasdynamic channel flow may be described by the following 


equations: 


energy, puC,T’ = up’ + (j2/c) (1) 
momentum, puu’ = —p’+  jB (2) 
continuity, puA = C (3) 
State, pb = pRT (4) 


Here, p is the density, u is the x directed velocity, C, is the specific 
heat at constant pressure, 7 is the temperature, 7 is the current 
density which is normal to the flow direction, and B is the mag- 
netic field which is normal to both current and flow direction 
o is the conductivity, which is assumed scalar, so that the current 
and effective electric field are connected by Ohm’s law: 

j = o( E — uB) (5) 


The magnetic Reynolds Number is assumed sufficiently small 
that B may be regarded as prescribed. The applied electric field, 
the temperature, and the conductivity will be assumed constant 
Then, combining Eqs. (1) and (2) and using Eq. (4) to eliminate 
B, we find 


pb’ = —(j?/uc) = RTE(j/uu')’ (6) 


* Senior Research Fellow in Jet Propulsion. 
** Professor of Jet Propulsion and Mechanical Engineering 
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Let ¢ j/u2u’, then 
¢’/e? = —[u%(u')?/oRTE] 


and, integrating from the origin, 
*Xx 
(1/¢) — [1/¢(0)] = (1/eRTE) | u*(u’ )2dx Gj 
J 
If the behavior of u is specified, that of 7 is readily calculated from 
Eq. (7), and that of B from Eq. (5). 
Particularly simple results are obtained when wu varies as a 
power of x. Suppose that 
“= ax” (8) 


I 


then, since the temperature is constant, the Mach Number, ./, 
ax” (4 RT)}/2 


, 


is given by 7 From Eq. (7 


1/¢ = (a'n?/oRTE) [x5"—'/(5n — 1) (9 


and, from the definition of ¢, 


j (¢RTE/a?n)(5n — 1)x~2" oE(5n — 1)/nyM? (10 


The pressure may be computed from Eggs. (1), (2), and (4): 

p |o( RTE)?/a5| [(5n — 1)/n2?}]x'~5" (11 
and the flow area from Eq. (3): 

A = (Ca‘n?/cE?RT)|[x™~'/(5n — 1) (12) 


Cis the constant value of puA. Finally, from Eq. (5), 


B = (E/u)\1 — [(5n — 1)/nyM?]} (13) 


In the special case of n = 1/5, these results are modified slightly, 
since then the integral in Eq. (7) is logarithmic 

All positive values of m greater than 1/5 yield flows in which 
continuous acceleration from zero velocity to an arbitrarily high 
velocity is possible, in a channel of continuously changing area 
The area and velocity variations are such that the conditions 
given by Resler and Sears! for smooth passage through the sonic 
For n = 1/4, the 
area is constant, and for this case, the solution passes through the 
singular point at M = 1, u = [(y — 1)/y](E/B). For n less 


velocity are satisfied for all such values of n 


than 1/4, the area decreases with increasing velocity, and for n 
greater than 1/4, it increases with increasing velocity. In view 
of the very rapid variation of both pressure and area with x for n 
very different from 1/4, values of m near 1/4 are likely to be the 
most interesting for plasma accelerators. 
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Experimental Evaluation of Heat Transfer 
With Transpiration Cooling in a Turbulent 
Boundary Layer at M = 3.2 


E. Roy Bartle* and Bernard M. Leadon** 
Convair Scientific Research Laboratory, San Diego, Calif. 
July 30, 1959 


 agetencrsesnes turbulent boundary-layer transpiration-cooling 
experiments’ yielded results which tend to correlate well 
with Rubesin’s theory.2, However, those results include a 
large radiation correction of uncertain accuracy, and it was 
assumed that the flow was two-dimensional although the channel 

* Staff Scientist. 

** Senior Staff Scientist 

The authors are grateful to the staff of the University of Southern Cali 
fornia Engineering Center for their cooperation during the use of their wind 
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was very small. To improve upon that study, it was concluded 
that heat transfer by radiation and conduction should be elimi- 
nated and the injection distribution made controllable during 
the tests 

Che new model, subsequently designed, has produced measure- 
ments that are believed to be more accurate than any previously 


reported. The test surface was a 5'/2 X 16 in. porous plate of 
sintered 5-micro-in. stainless-steel powder imbedded in the 
sidewall of a 3!/2 X 6-in. supersonic wind tunnel (M = 3.2) 


The model was divided into eleven compartments to permit 
adjusting the injection distribution to that necessary for 
a uniform surface-temperature distribution. Guard heating 
and cooling, together with highly polished, low-emissivity walls, 
reduced conduction and radiation losses practically to zero 
(A picture and further discussion of the model may be found 
in reference 9.) The heat balance on a segment of the porous 
plate, with the assumption that the coolant temperature attains 
the wall temperature at its outer surface, is then given by the 
equation 

St( pu Cp)i(Tw — Tr) = (pv Cp)w( Tw — Te) (1) 
in which St, T, p, v and u, and C, signify Stanton Number, 
absolute temperature, density, velocity components, and specific 
heat at constant pressure, respectively. In addition, the sub- 
scripts w, c, 1, and r refer respectively to wall, stored-coolant, 
free-stream, and adiabatic conditions. 

The measurements were made in the following manner: the 
tunnel was started, the desired stagnation conditions were set, 
and precooled or preheated gaseous nitrogen was metered 
through the porous plate at a constant flow rate. After thermal 
equilibrium was established approximately, the guard heaters 
and/or coolers were adjusted to eliminate conduction and 
radiation effects, and the injection distribution was adjusted 
to give uniform, steady, wall and coolant temperatures over 
the middle nine compartments, 14 in. Test-surface tempera- 
tures and coolant temperatures at a position 1/8 in. below the 
lower surface were then recorded. 

For each blowing rate, the coolant temperature was adjusted 
to three successive values near the recovery temperature and, 
in addition, to one at a nearly constant wall temperature of 
140°F. A plot of T,. — T. vs. Ty then gives 7, = 7, when 
T» — T. = 0. The recovery factor, o, is then given by (7; — 
T,)/(T° — T,), where 7,° is the free-stream total temperature 

The zero-blowing recovery factor, a, determined in a pre- 
liminary run using a solid plate, agrees well with the estimate 
o, = (Pr)'!/3, Also in the zero-blowing condition, but at 7, = 
140°F., values of the skin-friction coefficient were obtained 
by using a surface impact tube and the calibration curves of 
Fenter and Stalmach.5 The modified Reynolds analogy® then 
yielded values of St,, which agreed within 5 per cent with those 
of other investigators.4 

All heat-transfer data were obtained by running the model in 
a steady-state condition at a nearly constant wall temperature 


























Fic. 1. Heat-transfer reduction with nitrogen injection 


FORUM 79 















































LOs 
THEORY REF. 2 
| +s 
9 
FAIRED CURVE | | 
8 | 
£ M, = 3.20 
% 7 
: © —Reg = 13,100 
0 —Re, = 13,900 
6 + 
5 


Oo ! 2 3 4 5 6 7 
(PV, 1 
(Pu), St, 


Fic. 2. Measured recovery factor compared with theory. 






+ a 5 a 4 - 


— THEORY REF. 2 
/-—FAIRED CURVE 
a ; 


+ 





ef 
a 
“0 2 3 4 5 6 
(PV)y. 
(OU), Sty 


Fic. 3. Measured Stanton Numbers compared with theory 


of 140°F., and using Eq. (1). The local Stanton Numbers were ad- 
justed to 7/7, = 3.3 by a slight correction using a linear extra- 
polation of Deissler and Loeffler’s‘ results; they are presented in 
Fig. 1. Only data taken at compartments 8 and 9 are included, to 
ensure being far enough downstream that residual effects due to 
the transition from the solid nozzle wall (no injection) to the 
imbedded porous test surface (injection) had disappeared 


RESULTS AND DISCUSSION 


The blowing parameter employed is (pv)./(pu)St., which 
minimizes estimated? effects of Mach Number, Reynolds 
Number, and wall to free-steam total temperature ratio. The 
ratio o/a,, plotted in Fig. 2, is shown to decrease monotonically 
with increase of blowing in contrast to Rubesin’s estimate,? 
which he questioned at the time. 

The present local Stanton Number ratios, Fig. 3, fall above 
Rubesin’s? theory, indicating that it is an optimistic estimate of 
the effect of transpiration upon heat transfer. Present results 
also indicate that the earlier results of Leadon and Scott were 
perhaps overcorrected for radiation by as much as 25 per cent at 
the higher blowing rates. The data of Rubesin, ef a/.,° taken in a 
preliminary experiment, scattered seriously; but the average 
heat-transfer ratios form a smooth curve which falls slightly 
above present results at moderate blowing rates to as much as 
20 per cent above at the higher blowing rates. 

Recent skin-friction data reported by Pappas® fall above the 
present heat-transfer data by as much as 30 per cent at moderate 
to high blowing rates, indicating a strong effect of blowing upon 
the Reynolds analogy factor and one which is opposite in sign 
to that predicted by Rubesin.* It would be incautious to draw 
conclusions from this comparison, however, since part of the 
discrepancy may be due to dissimilar blowing distributions which 
could introduce the need for a cone-to-plate transformation factor 
applicable to the blowing parameter. 


These results are only a portion of the test program. A com 
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plete report covering the effect of injection on velocity profiles 
and Mach Number variations will follow. 
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On Three-Dimensional Boundary Layers 
Having Similar Solutions; 


M. G. Hall* 
Royal Aircraft Establishment, Farnborough, Hants, England 
August 29, 1959 


i and Hansen and Herzig?~* have investigated three- 
dimensional laminar incompressible boundary layers hav- 


ing similar solutions. Such boundary layers will, of course, have 
properties differing from those of two-dimensional boundary 
layers having similar solutions. The purpose of this note is to 
point out certain properties which, perhaps, would not be ex- 
pected simply from the addition of the extra dimension. For 
simplicity, we confine ourselves to similarity with respect to 
rectangular Cartesian coordinates (x, y, z) with the body surface 
represented by the plane z = 0, but the extension to orthogonal 
curvilinear coordinates is straightforward. The velocity com- 
ponents in the boundary layer are denoted by u, v, and w, with 
u = Uandv = V at the edge of the layer. The fundamental 
condition for similar solutions may be written 
u = U(x, y)f(r) 
v = V(x, y)g(¢) 
where ¢, the similarity parameter, is a scaled measure of the dis- 
tance along the normal to z = 0. 
We note first that the slope of a streamline within the bound- 

ary layer is given by 

dy/dx = v/u = (V/U)(g/f 
Thus, apart from the scale factor V/U and the scale factor in ¢, 
this slope will not vary from station to station. It follows that 

lim v/u = (V/U) X const. 

<0 
Therefore, since the streamlines at the outer edge of the bound- 


Tt The author is indebted to Dr. R. E. Meyer for helpful suggestions and 
the University of Sydney for the research grant which enabled this work to 
be carried out. 
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ary layer cannot form an envelope, neither can the limiting 
streamlines on the body surface; and we have the important 
result that the similar solutions do not admit boundary layers 
that separate locally along an envelope of limiting streamlines. ** 
It also follows that the total change in slope, 
tan-'[(V/U)}lim (g/f) — 1} + {1 + (V/U)? lim (g/f) 
t—0 <0 
will vary from station to station in more than scale alone. But 


the rate of change of slope 
(d/d¢)((v/u) + (V/U)] = (d/de\(g/f) 


will vary in scale only. 
The gradient of the velocity magnitude g = (u? + v?)!/? along 


a surface normal is given by 
(d/d¢)(q/k) = (U/k) ff’ + (V 


where k(x, y) is a scale factor. Clearly, no choice of the scale 
factor can in general produce a scaled velocity gradient that is 
independent of x and y._ For a two-dimensional boundary layer, 
however, we may put k = U, and the derivative reduces to 


(d/dg)(u/U) = f’ 


which is independent of x. Thus, in contrast to the two-dimen- 
sional case, the velocity gradients at different stations in a three- 
dimensional boundary layer having similar solutions differ in 
more than scale. 

The ratios to each other of the inertia, pressure, and viscous 
forces within the boundary layer can be obtained from the 
boundary layer equations in f and g. It can be shown that, for 
a fixed ¢, the ratios vary in general from station to station but 
will be the same for all stations having the same streamline slope 
V/U. Note that all stations in a two-dimensional boundary 
layer having similar solutions are in the latter category. 

An additional insight may be gained by rotating the coordi- 
nate axes in a given boundary layer and examining the profiles 
of the new velocity components #, 5. The new profiles in gen- 
eral will not be similar. However, if 7 = 0 locally, then it is 
found that 


v/(Vceosa) = g —f 


where a is the angle rotated through. Because V cos a may be 
taken as a scale factor, all profiles of such a component 0, which 
may be called the cross flow, will be similar. Evidently, this 
result would still hold if the x and y axes were only locally orthog- 
onal; if a three-dimensional boundary layer has similar solu- 
tions with respect to any orthogonal system of curves on the 
surface, the cross flow profiles will be similar. 
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** It can readily be shown that such local separation is impossible also 
with boundary layers having the solution for v compounded of similar solu- 
tions. An example is provided by the boundary layers of Hansen and Herzig’ 
under streamlines that are translates. 
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